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A detailed discussion of the calculation of rate constants for hydrogen atom transfer reactions based on the 
BEBO method is presented. Linear transition state models are used. A computer program using this method for 
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1. Introduction 

The bond-energy-bond-order (BEBO) method is a procedure for calculating the activation energies of 
hydrogen transfer reactions from bond energies. When combined with absolute rate theory, it also yields 
values for the rate constants. It was formulated over 10 years ago by Johnston and Parr [l], 1 and has since 
been applied with considerable success to the calculation of a large number of activation energies. Less fre- 
quently, it has been used to evaluate rate constants. Although the details of the BEBO method itself have 
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been published by Johnston [2], this aspect represents only a relatively small part of a rate constant calcula- 
tion. The purpose of this report is to give a detailed account, not only of the BEBO method and its 
theoretical background, but also of the absolute rate theory portion of the calculation. In addition, instruc- 
tions are provided for the use of a computer program which calculates rate constants based on the BEBO 
method. The discussion is limited to linear transition state models. 

2. Theory 

2, 1. Absolute Rate Theory & Transition State Model for BEBO Calculations 

For a bimolecular reaction, A + B — [AB]* — products, absolute rate theory utilizes the concept of a 
molecular complex made up of the two reactants. This complex is assumed to be in equilibrium with these 
reactants. The resulting expression for the classical rate constant k., is 

_ kT Q\,e-v»r 
Kt ~~ h & Q? t Ki} 

where k is the Boltzmann constant, Tis the absolute temperature, h is Planck's constant, Q? t and Qfi are the 
classical partition functions per unit volume for reactants^ and B, Q* is the classical partition function per 
unit volume for the complex, and V* is the potential energy of the complex relative to that of the reactants. 
The complex contains one unstable vibrational mode whose evolution brings about its dissociation into 
product fragments. The partition function Q* is evaluated with this mode missing. A detailed derivation of 
eq (1) which explains all its inherent assumptions has been given by Mahan [3]. Quantum mechanical correc- 
tions to the partition functions at room temperature and above need be applied only to vibrational factors. 
For a particular vibration of frequency v { , the quantum correction T,- is given by the expression 

Y { = — __L — ^^ where u,- = hvJkT (2) 

sinh (u,/2) 

We assume that all vibrational modes are independent so that the total quantum correction for a particular 
species is simply the product of terms given by eq (2), one for each vibrational mode. There is also a quan- 
tum correction to the unstable vibrational mode of the complex which we denote by Y*. This results from 
the effect of quantum mechanical tunneling through the potential barrier between reactants and products. 
It will be considered in detail in section 2.5. Applying these quantum corrections to eq (1) yields the rate 
expression 

kT q: {?rnv*e- v '>« 

h Qt, {prn Qf, {prf} w 

The general class of reactions we are considering has the form 

A-H + £• - A~H~B -A* + H-B (4) 

Radical B* abstracts a hydrogen atom attached to A, the net result being the transfer of H from A to B. For 
this system, we take the most general transition state to be linear, having up to 5 mass points. Its structure 
and the notation which we shall use are shown in figure la. There can be up to four internuclear distances, 
/?,, R b , R„ and R d . The bonds associated with R a and R d will be assumed to be rigid. (The two vibrational 
modes involving these bonds will have infinite frequencies and need not be formally included in the calcula- 
tions.) Thus, there are only two vibrational stretching modes to be considered for this molecule, one of which 
will be unstable. These modes arise from the stretching of the two central bonds b and c which are shown by 
dotted lines to indicate their unstable character. Of the five possible masses, M 3 will normally be that of the 
hydrogen atom; the other masses will be assigned values in the manner described below. The three angles 
^j, ^f 3 , and ¥ 4 are defined by the bonds (a,6),(o,c), and (c,d) in the plane of the figure while the primed sym- 
bols denote the corresponding angles in the plane perpendicular to the figure. Changes in these angles from 
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R a Rb Rc Rd 

M| — M2«"'M3' ,, 'M 4 — M 5 

% % % 

FIGURE la. Notation for five mass point linear transition state. 

180° give rise to three doubly degenerate bending vibrations. To calculate the frequencies needed in eq (3), 
we require values for the two stretching force constants associated with bonds b and c, and three bending 
force constants arising from the three bond angles. As we shall see, these values can be generated by the 
BEBO process. 

Within the framework of the transition state structure shown in figure la, it is possible to include all types 
of reactions implied by eq (4) by considering four cases; one having a 3 point transition state, two having 4 
point states, and one having a 5 point state. These four cases are shown in figure lb. In this figure, the 
subscript s appearing on the internuclear distances and force constants denote equilibrium values found in 
reactants or products. Because bonds a and d are assumed to be rigid, their bond distances will always be 
denoted by the single symbols /?„ and R&, respectively. The bond distance between M s and M 3 goes from 
R b , to R b in the transition state, while that between M 3 and M 4 goes from o° to R c in the transition state. In 
the transition state, the force constant F bs is modified and combined with that of the newly formed bond be- 
tween M 3 and M 4 to produce two force constants F„ and F a . F„ corresponds to the stable symmetric stretch 
and F„ to the unstable asymmetric stretch. In cases IVa and V, the bending force constant F*2i becomes F+j 
in the complex. The newly formed bond angle made by M z , M 3i and M 4 leads to the force constant F* 3 in all 
cases. Finally, in cases IVb and V, we also have an additional bending force constant F+ 4 which goes to F* 4 , 
in the second product. The force constants associated with the out-of-plane bends are not shown since they 
are the same as the in-plane constants. 

The way I have chosen to assign values to the mass points is somewhat arbitrary and is best explained by 
an example. Consider the reaction 

CH 3 -CH 2 -H + CH 3 - CH 3 -CH 2 ..H..CH 3 - CH 3 -CH 2 . + H-CH 3 
Species A B CD 

which is the abstraction of hydrogen from ethane by methyl radicals. The masses are assigned according to 
the following rules: 

1) The mass of the transferred H is always assigned to M 3 ; therefore M 3 = 1.008 atomic mass units 
(a.m.u.). 

2) The mass of the atom joined to the transferred H in reactant A is assigned to Af 2 ; in this case M z — 
12.011 a.m.u. 

3) The masses of all the remaining atoms in A are added and assigned to Afjj thus in this example M x = 
17.051 a.m.u. 

4) The mass of the atom joined to the transferred H in the product D is assigned to M 4 ; here M 4 = 12.011 
a.m.u. 

5) The masses of all the remaining atoms in D are added and assigned to M s ; thus M s = 3.024 a.m.u. in 
this example. 

Different models for the transition state, and different ways of arranging the masses in linear models have 
been explored in a limited way by Johnston [4] and by Sharp & Johnston [5], They did find significant dif- 
ferences between various options. Presumably, complete vibrational analyses of the reactant and complex 
would yield more accurate rate constants than the linear models outlined above. Unfortunately, complete 
analyses are extremely complex even for fairly small molecules, and the ability to program the calculations 
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FlCURE lb. Reaction cases to be used in BEBO calculations. All transition states shown here are linear. Masses are denoted by M, inter- 
nuclear distances by R, and force constants by F. The subscript s denotes bond distances and force constants in the stable reactants and 
products. 

in a general manner would be lost by such an approach. Also, it is unlikely that all of the force constant 
values required would be available for a complete analysis. In view of the crudity of the rest of the calcula- 
tion, it is unnecessary to strive for high accuracy in the vibrational analysis. Intuitively, one expects that the 
major features of these reactions are controlled by the nature of the atoms adjacent to the H atom being 
transferred, with the effects from the remainder of the molecule appearing in the bond energy values. If this 
is the case, then the linear models should at least be able to match trends within homologous series. 

So far, we have seen in this section that evaluation of rate constants by the use of eq (3), based on the 
linear models shown in figure lb, requires a knowledge of the potential energy V* of the complex, two 
stretching force constants, and from one to three bending force constants. The potential energy of all of 
these linear models could, if it were known, be shown on a 2-dimensional contour diagram like that shown in 
figure 2 where the independent variables are the bond distances R k and i? c . The required value of the poten- 
tial energy Vl$ that at the saddle point position shown by the asterisk. For a region close to the saddle point, 

608 




FIGURE 2. Typical potential energy diagram for H atom exchange reaction. 
The position of the saddle point is shown by the asterisk The direction q is 
that in which the potential energy decreases most rapidly. The direction o 
is perpendicular to the q direction. 



it is customary to assume that the first derivatives of V with respect to R b and R c are negligible, and that the 
potential energy can be approximated by a power series containing only quadratic terms. Thus, for small 
displacements from the saddle point, we have 



where F„ = 



B 2 V 
dRV 



F = F = 



25 V * F n (8R b ) 
d 2 V 



dR b dR c 



F = 



+ 2F 12 (8R b )(5R c ) + F^bRy 
d 2 V 



(5a) 



dR 2 



These derivatives are evaluated at the saddle point, and are, by definition, the stretching force constants of 
the complex. In matrix notation, this equation is 



25 V « (5R)fF r (5R) 
where F r = [5» 5» ] and 5R = [ *£ ] . 



(5b) 



This is the force constant matrix that will be used to calculate the vibrational stretching frequencies. 

Starting at the saddle point, suppose we move in the direction in which V decreases most rapidly; call this 
the q direction, and let a denote the direction perpendicular to q. These directions define a rotated set of 
cartesian coordinates which we assume makes an angle a with the R b axis; (positive a is measured in the 
counter-clockwide direction). The transformation between the two sets of coordinates is given by the equa- 
tion 



R 



L r J L 



Rc 



cosa 
sina 



cosa 



]m 



UP 



(6) 



the matrix U can now be used to express changes in Fat the saddle point in terms of changes in q and a in- 
stead of R b and /? c . Thus, eq (5b) becomes 



25 V * (5R)fF,(5R) = (U5P)fF r (U5P) = (5P)«UtF r U)(5P) 



(7) 
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The matrix UfF r U has the elements 



(UtF,U) n = F n cos 2 a + 2F 12 cosasma + F 22 sin z a 
(UtF,U) 12 = (UtF,U) 21 = (F 22 - F n )cosasina + F 12 (cos 2 a - sin z a) 
(UtF r U) 22 = F u sin 2 a — 2F 12 cosasina + F 22 cos 2 a: 



(8) 



As we shall see in the next section, the BEBO method provides values for the second derivatives of F(i.e., 
the force constants) in the q and a directions. This will allow us to evaluate the matrix UfF r U. The stretch- 
ing force constant matrix F r , can then be obtained by inverting the transformation given by eq (6). 

In this section I have presented a formula (eq (3)) for the rate constant and outlined the factors required to 
evaluate it. The details of the BEBO method will be given next. It will provide values for V* and all of the 
necessary force constants, both the stretching and the bending ones. 

2. BEBO Method 

The BEBO method is based on the concept of bond order. In the reactants the bond b of figure la is said 
to have a bond order of unity, while in the products, its bond order is zero. The reverse of this situation 
holds for bond c. BEBO assumes that during the reaction, the total bond order of the two bonds is con- 
served; if n is the order of bond b, and m of bond c, then we have always n + m — 1. This is the basic 
assumption of the method. One bond is breaking at the same time that the other is forming. To apply this 
conservation condition it is necessary to relate the energies and lengths of bonds b and c to their bond 
orders n and m. 

For the relationship between order and length, Pauling [6] proposed the formula 



R n = R, - \ln{n) 



(9) 



where R, is the length of the bond which is considered to be representative of a single bond between the two 
elements of interest. The parameter X is taken to have the same value for all element pairs. A plot of bond 
length versus the logarithm of the bond order is shown in figure 3 for certain element pairs. The data were 
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obtained from table 4-3 of Johnston's book [2]. Values of X obtained from least squares fits are given in the 
figure for the different bonds. The constancy of X does not seem particularly striking. Pauling chose a value 
of 0.26 for X; he writes, "This equation, which is based upon the study of interatomic distances for non- 
resonating and resonating covalent bonds in simple non-metallic substances of known structure, is found to 
agree reasonably well with those data for metallic crystals which are suited to a check on its validity, and its 
use permits a penetrating analysis of the structure of metals and intermetallic compounds to be made. There 
is some evidence that the constant . . . varies with the kind of atom and with the type of bond; but the 
evidence is not sufficiently extensive to lead to the determination of the nature of this variation." Certainly 
Pauling's value doesn't appear to have been based very heavily on the data in figure 3 since none of these X 
values are close to 0.26. Although 0.26 can hardly be construed as universal, it has nevertheless been the 
value used for most BEBO calculations. There appears to be no reason why a different value shouldn't be 
used if it gave better results. 

Consider next the dependence of bond energy on bond order. Johnston [4] proposed the following rela- 
tionship between the two quantities 



E n = Eji> 



(10) 



where E, is the bond energy of a single bond and is analogous to R s of eq (9). Note that this energy is the 
electronic dissociation energy of the bond in question; the zero point energy is not meant to be included in 
E 5 . Plots of ln(E) versus \n(n) are shown in figure 4 for the same bonds used in figure 3. The data are again 
from table 4-3 of Johnston [2]. We see that p depends on the kind of atoms in the bond. If more than one 
bond type occurs for a pair of atoms, then it is possible to extract values for p from plots like figure 4 pro- 
vided we are not unduly bothered by a lack of linearity. When only a single bond type exists, then some 
other method must be devised. Actually, since we are interested in E n and R n for bond orders less than unity, 
even if multiple bonds were available for a plot like figure 4, some method of extrapolating to zero n would 
be necessary. Johnston [2], inspired by Badger's rule for force constants, has devised a way. Let us first elim- 
inate n between eqs (9) and (10); this yields 



ln(EJE.) = (pl\XR, - Rn) 



(ID 



This expression is analogous to Badger's rule (see Herschbach & Laurie [7]), which is a universal empirical 




FIGURE 4. Plots of lr(£J versus ln(n)for certain bonds. 
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relation having the form r = a u - b tJ log(/), where r is the bond distance, f its force constant, a y and b u are 
constants, and i and j are the numbers of the rows in the periodic table in which the bound atoms are 
located. Johnston [8] found that plots of log(/) versus r extrapolated very nicely to two-atom Lennard-Jones 
noble gas clusters. For clusters having Lennard-Jones parameters a and elk, the "bond" distance is 2^tr and 
the "force" constant is 4G.06(e/&)/a 2 . He then examined plots of ln(E„) versus R„ to see if a comparable 
extrapolation would be possible. The results are shown in figure 5. The data are mostly from Johnston [2], 
tables 4-3 and 4-1. Values of E„ and R„ for the He-Ne cluster were taken from Gilliom [9]. The energies for 
the bonds examined in figures 3 and 4 are supposed to extrapolate to the Ne-Ne cluster. The lines shown 
were drawn to connect the corresponding single bonds with this cluster. Points corresponding to multiple 
bonds fall more or less in the general direction of these lines. The assumption made in BEBO is that such an 
extrapolation adequately represents the bond energies for n < 1. Therefore, if we have a bond A-H, where A 
is some atom in the first row of the periodic table connected to an H atom, and R, and E, are its bond length 
and energy, then if this bond were perturbed in some fashion so that its bond length were greater than R„ 
then its bond energy would fall on the line drawn between the A-H and He-Ne points. Bonds involving atoms 
A from other rows of the periodic table will extrapolate to the appropriate rare gas-helium cluster. The slope 
of the line joining A-H to the cluster is, from eq (11), -p/X. Since the value of X has been chosen, we have a 
way of calculating p for the A-H bond of interest. Formally, in this case, 



p = 



R, — Rb 



ln{E H ^ Nt lF,) 



(12) 



The parameter p thus depends onX, the bond energy and internuclear distance of A-H, and the interaction 
parameters for the appropriate rare gas cluster. 




FIGURE 5. Extrapolation of bond energy to large bond distances. A of A-H 
is an atom in the first row of the periodic table in this erne. 

We have now almost all of the information needed for the BEBO calculations. Consider a triatomic com- 
plex A-H-B; there are three interactions; two between H and the atoms A and B considered above, and the 
interaction between A and B themselves. If H is to form stable bonds its electron spin must be opposite each 
of the spins of A and B. Consequently, A and B will have parallel spins and must repel each other. Johnston 
uses one half the value of the Sato [10] triplet function to represent this interaction. He uses the modified 
function because it more closely approximates the calculated H-H triplet interaction. This function has the 
form 



612 



V, = E u E{\ + E) (13) 

where E = V6e"* A, S AR, = R t - R a = R b + R e - R u . E„ is the electronic dissociation energy, R„ the 
equilibrium internuclear distance, and j8 the Morse parameter (see Herzberg[ll] p. 101) of the ground state 
of the diatomic molecule made up of A and B. Values of these paramaters for a number of such atom pairs 
are given in Table 1 1-1 of Johnston [2], AR, is the difference between the actual distance R, between A and 
B in the complex and the equilibrium distance R„ it would have in the diatomic molecule. It is worth point- 
ing out that many people use E„ as an adjustable parameter to fit the BEBO calculations to their experimen- 
tal data. Other forms of the triplet function have been used and are discussed briefly in the Appendix. V, can 
be expressed as a function of n, the bond order of the 6 bond, through the conservation condition n + m — 
1, and through eq (9) which gives the distances R b and R c in terms of n and m. 

We are now able to give the BEBO expression for the energy of the complex in terms of the bond order n. 
The energy is assumed to be given by 

V[n) = E b , - E b ,n" - £„m' + V t (n) = EJ} - n") - E a {\ - n)" + V t {n) (14) 

E u and E a are the single bond energies (electronic) for bonds b and c, and the parameters p and q are calcu* 
lated from eq (12) for 6 and c, respectively. When n — 1, then m— -0, V t — 0, and F— 0, so that the energy is 
measured relative to the energy of the reactants. When n— -0, then m—1, V t — 0, and V~*E b ,~E a which is 
the difference in the bond energies. BEBO assumes that the maximum value of Fin the range 1 2: n > is 
the desired potential energy of the saddle point. This value V*, is obtained by substituting into eq (14) that 
value of n which makes dV/dn = 0. In what follows, all quantities are considered to be evaluated at the sad- 
dle point. 

Next, we must determine the stretching force constants in the q and a directions shown in figure 3. Equa- 
tion (14) does not give the complete potential surface, but only that portion lying along the line of constant 
total bond order. BEBO assumes that at the saddle point, this path of constant bond order lies in the q 
direction. This assumption will enable us to calculate the force constant F g a d*V/dQ 2 from the second 
derivative of Fwith respect to n, which we get by differentiating eq (14). 

From eq (9), we can calculate the changes produced in R b and R e when n is changed. In vector notation 
these are 

SR = [^] = - X [ -i7»] 6n;m = 1 - n - (15) 

Because a change in n for constant total bond order is supposed to produce a move in the Q direction, the 
slope of a line in this direction can be gotten from eq (15). It is 

5R c 15R b = -nlm = tana (16) 

where a is the angle which q makes with the R b axis as discussed earlier. From eq (15) we can show that 

cosa = m/^J{n s + m 2 ); sinor = -n/yj{n 2 + m 2 ) (17) 

The matrix U defined in eq (6) can now be written in terms of n and m. 

U = V(« 2 + m 2 ) I ~ n m \ (18) 

By means of eq (6), 5R can be expressed in terms of 5P; i.e., 6q and 5a. Combining the differential form of 
eq (6) with eq (15), we get 

U6P=-\[_ 1 (; m ]5a (19) 
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Solving for 6P gives 






(20) 



As expected, a does not change when n changes. From eq (20), we have for the derivative of n with respect to 
Q 

dn dn 1 nm 



bg dg 



X V(n 2 + m 2 ) 



(21) 



The second derivative of V with respect to g is obtained from the sequence 

dV dV dn 

dg dn dg 



<PV <PV 



dg 2 



m 



2 dV d*n 
dn 2 \ dg / dn dg 2 



Since dV/dn = at the saddle point, we have 

d 2 V d 2 V n 2 m 2 



F a =- 



dg 2 dn 2 Wn' + m 2 ) 



(22) 



This gives one of the stretching force constants. 

In the a direction, the stretching motion is assumed to be that of a normal molecule. Thus Badger's rule 
should be applicable. This says that the bond distance is proportional to the logarithm of the force constant, 
while eq (9) says that the bond distance is proportional to the logarithm of the bond order. Therefore, the 
force constant should be proportional to the bond order. We assume that 



F b =F bl n, and F e =F a m 



(23) 



where F bs and F a are single bond force constants. Consider the change in fwhen R b and R c are changed by 
motion in the a direction. This is assumed to be given by 



2(dV) a =F bl n(bR b )l+F cl m(5R c )l+^- (bR,)l=Fj[6oY 



(24) 



To evaluate F„, we must express (6R b )^ (5i? c )5, and (8RX in terms of (5a) 2 . From Eqs (6) and (18), we have 



For 5^=0, 



Therefore, 



6R = 



8R b 
8R C 



1 



y/(n 2 +m 2 ) 



m n 



— n m 



(8R b ) a = nbal V(n 2 + m 2 ), {6R c ) a = mbalyj{n 2 + m 2 ) 
(bRX = (bR b ) a + (bR c ) a =(n + m)balyj{n 2 + m 2 ) = boly](n 2 + m 2 ). 



2{bV) =(F bt n>+F a m 3 +^f) {boYI(n 2 + m 2 ). 



(25) 
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Comparing this with eq (24) gives 



Fc = 



F bs n 3 + F ts m 3 + PVJdR 2 
n 2 + m 2 



(26) 



The method assumes that if V is expanded at the saddle point in terms of q and a then there is no cross 
term; i.e., d 2 Vldgdo is assumed to be zero. Thus, we have 



The use of eq (7) shows that 



Inverting this equation gives 



25 V = (5P)f 



UfF r U = 



F e 
F a 



F e 

Fa 



(5P) 



F r = (Uf)- 



F e 

Fa 



u-» =u 



F e 

Fa 



ut 



(27) 



(28) 



(29) 



where use has been made of the fact that U -1 = Uf. Substituting eq (18) into (29) gives the desired stretch- 
ing force constant matrix. 



F r = 



1 



n 2 + m 2 



m n 
— n m 




F e 

Fa 




m —n 
n m 



n 2 + m 2 



F e m 2 + F a n 2 , —F e mn + Foinn 
—F Q mn + F a mn, F e n 2 + F a m 2 



^11 ^*12 
^21 ^22 



(30) 



To complete the discussion of the BEBO method the bending force constants will now be evaluated. Con- 
sider first the one involving M 3 as the center mass. This will be F+ and appears in all of the transition states 
shown in figure 2. It is defined as the second partial derivative of V with respect to the angle made by the 
bonds 6 and c, with the bond lengths Rf, and R c held fixed. At equilibrium, this angle is 180° for our transi- 
tion state models. The geometry, when the angle is less than 180° is shown in figure 6. To get F+ , we dif- 
ferentiate V twice, 



(— ) 

( d 2 V \ 

V W ! R b ,R c 



dV, dV, dR t 



d<t> 

d 2 V t 
d<}> 2 



dR t d4> 

d 2 V t (dR t Y 
dR 2 V d<t> } 



dV, d 2 R, 
dR t d<}> 2 



(31) 




FIGURE 6. Definition of center bond angle. 

The derivatives of V, with respect to R, can be gotten from eq (13). The dependence of R, on can be deter- 
mined from the following vector relationships, 
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Rj • R c = i? 6 i? c cos0 
R t — R c — R* 

/?? = R, • R, = Rl + Rj - 2R b Rc cos0 
dR, _ R b R c s . n0 _ Q for = 18Q0 



30 R, 

d 2 R t _ R b R c cos0 

d0 2 R, /?? 



75/5 sin 2 0- -^for0 = 180* 



Thus, 



^Va = ~ 



ar f R b R c 



R, 



dV, R b R e 



(32) 



(33) 



dR t R t 3R t R b + R e 

The other two bending force constants F^ and iv 3 are assumed to obey Badger's rule. We assume 

^V 2 = F**Ji> and F* 4 = F* 4j m (34) 

This concludes the BEBO part of the calculation. It has provided us with the potential energy V* of the sad- 
dle point, the stretching force constants F 1V F^, and F 12 and the bending force constants F+ 2 , F+ 3 , and F^. 
In the next section we shall use these force constants to carry out a frequency analysis for each of the transi- 
tion state models shown in figure 2. 

2.3. Vibrational Analysis 

As we have seen in the force constant derivations, the potential energy V of the most general 5 mass point 
complex can be considered to depend on the variables R a , R b , R c , R d , ¥ 2 , ¥ 3 , ¥ 4 , ¥ 2 ', ¥ 3 , and ¥ 4 '. These are 
called the internal coordinates. Because our model is linear, V increases when any of the angles departs 
from 180°. Since we assume a and d to be rigid, R„ and R d need not be included in the list of variables. For 
the time being, however, they will be included in the analysis. Let F be the complete force constant matrix 
for the complex. We have 



F = 



Ml ^12 
**21 **22 







'*2 







*3 



#4 



K 



F4 



FJ. 



(35) 



The two infinite force constants come from the use of rigid bonds for a and d. Let S be the (column) vector 
which denotes small changes in the saddle point values of the variables. 



St = [5/?„, hR b , 6R C , 6R d , 6%, 8%, S¥ 4 , 8% 6% S¥ 4 '] 
The potential energy is assumed to be given by 

V - V* = i/ 2 StFS 



(36) 



(37) 
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Suppose there exists a matrix G, such that the kinetic energy in terms of the internal coordinates is 

T = V6StG-*S (38) 

Consider a new set of coordinates Q s the so-called normal coordinates, related to S by the linear transforma- 
tion 

S = LQ (39) 

such that 

V - V* = VfcQtAQ = AV (40) 

T = i/ 2 QtEQ (41) 

where A is a diagonal matrix having elements X„ and E is the identity matrix. In this coordinate system 
there are no cross terms in V and T, 
Let Qi denote the i'th normal coordinate. The Lagrangian equations of motion for the system are 

where L = T -AV = V*[QfEQ - QtAQ] = V 2 [ZQf - E\jQ]] , (43) 

d k- = Q h (44) 



dL 



- -KQi. (45) 



Therefore 



The solutions of this equation are 



dQi 
Q. + \.Q. = 0. (46) 

Q, = Q?cos(\?t + e,). (47) 



Thus the \? = 2irv t are the frequencies of the vibrations of the Q t coordinates. These are called the normal 
mode vibrations. 

Solving eq (39) for Q, and substituting into (40) and (41) yields 

V _ V * = i/ 2 (L-S)tA(L- 1 S) = i/ 2 St(L- 1 )tA(L- , )S (48) 

T = y 2 (L- J S)fE(L-'S) = VfcStfL-'JtEtf/'JS (49) 

Comparison with eqs (37) and (38) yields 

F = (L-)tA(L-) (50) 

LfFL = A (51) 

G 1 = (L-')tE(L-) (52) 

LfG'L = E (53) 
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Next, solve eq (53) for Lf = L _1 G, substitute this into eq (51) and multiply by L on the left. This gives 

GFL = HL = LA (54) 

as the set of equations which determine the transformation L, Written out, eq (54) is 

Lj[H v - 5 y X t ]I ifc = (55) 

This equation has solutions if the determinant 

|H - EX* | = (56) 

This is the so-called secular equation which must be solved to get the X*, the eigenvalues of H and thus the 
normal frequency values. Before doing this, it is first necessary to evaluate the matrix G. 

Equation (38) gives the kinetic energy in terms of the internal coordinates. As such, it does not include the 
kinetic energy of the center of mass or the rotational energy. We need to express the kinetic energy in terms 
of cartesian coordinates, transform the result to internal coordinates, and subtract out the center of mass 
and rotational energy. This will yield G" 1 . Let us begin by expressing the internal coordinates in terms of 
cartesian coordinates. Assume that the molecule lies along the x axis. A particular mass point M,- will have 
coordinates (x ( , y t , z t ) where y, and z,- are small and describe the departures of the molecule from linearity 
during bending vibrations. Because y t and z t are small, the bond distances can be expressed as functions of 
the %i only. Thus, 



Rb - * 3 - x 2 



R c = x 4 - x 3 



Rd — x* — 



(57) 



Since there are 5 cartesian x coordinates we need one more coordinate for the internal system. This is taken 
to be the ^-component of the center of mass of the molecule multiplied by the total mass, and is defined by 
the equation, 



Mx = E s Mix, 
where M — E 5 M t 



(58) 
(59) 



In matrix form these equations are 



R = [ Mx] = 



R. 




-1 


I 













*i 


R, 







-I 


1 










X. 


Rc 


= 








-1 


1 







x : 


R< 













-1 


I 




X 4 


Mx 




M x 


AT, 


M 3 


M t 


M s 




X 



= MX 



(60) 



Note that the vector R is basically that defined by eq (15). Here we have included R a and R d . 

We must next express the bond angles in terms of the cartesian coordinates. Consider ^ 2 , the angle 
formed by bonds a and L The geometry and notation for this angle are shown in figure 7. The two vectors 
along the bonds a and b are given by 
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rt 2i = [*i - Wi - y 2 ] 
rt 23 = fo - * 2 ,y 3 - y 2 ] 





*~ 


^ 






•r 


\£/ a> %/ 


^ 




x 


r \^ 




A 


3 



FlGVRE 7. Geometry of the bond angle V 2 



(70) 



¥ 2 is related to these by 

r 2i " r 23 = r 21 r 23 cos¥ 2 . (71) 

Substituting eq (70) into (71) gives 

-R a R b +(y 1 - y 2 )(y 3 - y 2 ) = {[R\ + {y x - y 2 )*][/?? + (y 3 - 72 )*]}*cos¥ 2 (72) 

Because the y t are small compared to R a and R b the radical can be expanded to give 

-R a R b +( Jl - y 2 ){y 3 - y 2 ) = [R a R b + i/ 2 (/? 6 //? )(y 1 - y 2 ) 2 + V 2 (RJR b )(y 3 - y 2 ) 2 ] 

• cos¥ 2 (73) 

Let ¥ 2 = 180° + 5¥ 2 where 5¥ 2 is small. Then 

cos¥ 2 = -cos(8%) « -1 + V 2 (8%f 
Substituting this into eq (73) and keeping terms through second order gives 

(5¥ 2 ) 2 = (URDfa - y 2 ) 2 + (2//?,,/?,)^ - y 2 ){y 3 - y 2 ) + (\IRl)(y 3 - y 2 Y 
5*2 = -[(/! - y 2 )IRa + iy 3 - y 2 yR„] = -yJR. + (UR. + UR>)y 2 - y 3 IR„ (74) 

To see why the minus sign is needed, let/! = y 2 = 0; then fory 3 > 0, ¥ 2 < 180°, so that 5^ 2 must be < 0. 
There are analogous equations for the angles ¥ 3 and ¥ 4 ; there is also a set, identical in form, for the angles 
¥/ in the x-z plane. These contain the z,- rather than the y, coordinates. In these equations, the equilibrium 
values of R a , . . . ,R d will be used. 

The set of equations typified by eq (74) gives 3 equations in terms of the 5 y ( coordinates; two more are 
needed. We have one defining the y coordinate of the center of mass, like eq (58), and another defining a 
quantity 77„ which is given by the equation 



77, = E'MiX'Yt 

77, is related to the z component of the angular momentum m, by the relation 

m z = 77, 

619 



(75) 



(76) 



The xf, are the equilibrium x t values; these can be gotten relative to the center of mass component x, by in- 
verting eq (60) and inserting equilibrium values for /?„,... ,R d . In matrix form, these equations relating y t to 
the bond angles in the x-y plane are, 



* - 







fyz 




V 




5^3 




y* 


= 


s<A 4 


= 


My 




*7» 






My 





~Qa Qa + Qt ~Qb 

~Q b Qb + Qc ~Qc 

—Q e Qc+Qd ~Qd 

Mjcl M^l M 4 xl M s x' 5 



M lX \ 



M 9 



M, 



M A 



M, 



y 3 
y* 

7s 



= AY (77) 



where g„, . . . ,Qd are the reciprocals of the equilibrium values of /?„, . . . ,R d . There is an analogous equation 
involving the z, coordinates. 

Having obtained expressions (60) and (77) for the internal coordinates in terms of the cartesian coor- 
dinates, we can nowjnvert these equations and insert them into the expression for the total kinetic energy 
which we shall call T. Therefore 



T = >/ 2 XfD m X + 1/2 YfD m Y + */ 2 ZtD m Z 

= V2 R-KM-^tD^M-'JR + 1 / 2 ^t(A- , )fD B ,(A- 1 )* + z-term 

= J/2RtG r -R + V 2 tiGi % t + z-term 

= l AT + 1 /2M(x z + y 2 + F) + Vkimj + mf)// 

= i/fcRtGr'R + VfetftG/ 1 * + z-term + \kM&* + y 2 + ?) + V 2 (m 2 y + rn^ll 
where / = E 5 M t xf is the moment of inertia, and 



(78) 



D m = 



M, 







Af, 







M 4 



M. 



(79) 



We can satisfy eq (78) by writing G r "' and G; 1 in the partitioned forms 

G r > = 



G r "' 

AT 



(80) 



G;» 











/-» 











AT 1 



GJ = 



We can get G r and G* simply by inverting G r -1 and G; 1 . 
This gives 



= G*: 1 



; ' - &' 





M 



Gj, — 



G, 
I 

M 



= MD m -»Mt 



= AD.; 1 Af 



(81) 
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Since D m is diagonal its inverse is easily evaluated and we therefore require only matrix multiplications to 
get G r and G+. 

The complete G matrix for the internal coordinates in partitioned form is 



G = 



F = 



G r 











\J^ 





J) 





G^» 


stant 


matrix, 


eq (35), 


~F r 





o — 





F* 











F,< 



(82) 



(83) 



Note that F r here is like eq (30), but contains the two infinite force constants corresponding to the rigid a 
and d bonds. The matrix H in partitioned form is 



H 



XJr* r 










H r 











G^F^ 





_ 





H^ 











G^F^' 










H^.» 



(84) 



Because H factors in this way, we can set up separate secular equations for the stretching and bending 
modes. Note that H is normally unsymmetric. 

Before solving the secular equations, let us write down explicit expressions for G r and G*. The direct 
evaluation of G r from eq (81) yields 



G r = 



l+/*2 


~H 











~/*a 


th+^3 


~lh 











~H 


Ih + P* 


-/*• 











-A* 


M 4 +Ms 

















M 



(85) 



where the fi { are the reciprocals of the masses M t . Comparison of this equation with eq (81) yields G r . 
Because we are treating the a and d bonds as rigid, the stretching part of the problems is equivalent to a 3 
mass point system where the first mass is M x + M 2 and the third is M 3 + Af 4 . The resulting 2x2 matrix is 
the one actually used in the calculation. It is 



G r (rigid end bonds) = 



Mi 



Ml+Mi 



-)* 



+ /*3 



~lh 



Ma 



/*3+M. 



.(- 



M* 



M4+MS 



(86) 



The stretching force constant matrix to be used with eq (86) is that F r as given by eq (30). 
The G^ matrix elements for this 5 point case are 

(G*)„ = Qlh + qIh 3 + (e« + Qb) 2 ^ 
(G#)„ = Qlh + e?M* + (e* + Qc)*ii a 

(G#) M = Q 2 cfh + Qd/i s + (Gc + Qd)*li t 

(G#) u = -e*[(c- + e*K + (q» + qM 

(G*) 23 = - Q c [(q„ + G«Vj + fe* + Qd)fi*]i (G*) 13 = QtQcfh 
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(87) 



There are the expressions used in the calculation. Actually, they were not derived from eq (81) but were ob- 
tained from Wilson et. al. [12]. However, eq (81) was used for a numerical check of eq (87). To get the matrix 
elements for the two 4 mass point cases, simply delete from eq (87) those elements which contain either a 
missing q or a missing /t or both. Do the same for the 3 point case, but delete also (G^) 13 ; (there is only one 
element, (G^) 22 , in this case). 

We are now ready to consider the secular equation. For the rate constant calculation only the X* are re- 
quired, so that a solution of eq (54) for the transformation matrix L is not necessary. Nevertheless, L is easily 
obtained and is convenient to have for the purpose of illustrating the actual vibrational motions of the com- 
plex. Thus we shall solve eq (54) as well as eq (56). According to eq (84), there are two secular equations to be 
solved (H^ and H^# are equal). Because we are using rigid a and d bonds, the dimension of H r is 2 X 2. The 
maximum dimension of H* is 3 X 3 and occurs for the 5 point model. Thus a solution of a 3 X 3 problem 
will suffice for our purpose and will also illustrate how an n X n problem is to be solved. 

We begin by assuming that eq (56) has been solved. In the present work this was accomplished by expand- 
ing (56) and solving the resulting polynomial in X. In our case, the maximum degree was cubic, so that this 
part of the calculation was easily performed. As eq (47) shows, the desired frequencies are v k = XJ /2 /27T. For 
the stretching modes of the complex one of the two frequencies will be imaginary because its X* value will be 
negative. As mentioned earlier, this corresponds to the asymmetric stretch. 

Consider now eq (55) for a general 3 X 3 H matrix. Written out in full, it is 



(H u — \ k )L lk + H 12 L 2k + H 13 L 3k = 
H 21 L lk + (i/^ — \ k )L 2k + H 23 L 3k = 
H 31 L lk + H 23 L 2k + (/7 33 - \ k )L 3k = 



(88) 



where X* is one of the three values of X determined from the solution of the cubic (in this case) eq (56). 
Divide the first two of these equations by L 3k , and define the ratios g ik = L ik IL 3k . This yields two equations 
to be solved for the two unknowns g lk and g 2k . 



(#n - X*)g u + H 12 g 2k = -H h 



(90) 



H 2 lglk + (#22 ~ \k)glk = -#23 

We get two g ik values for each value of X* substituted into eq (89), or six g ik values in all. Using these values, 
we can express L in terms of the product of two matrices defined by 



L = 



gll 


gll 


S\ 3 




L 31 








8n 


#22 


823 







L 32 





1 


1 


1 










L 



= r? 



To determine the components of t, insert eq (90) into eq (51). We get 

ftrtFTf=A=ftfrtFr 



(90) 



(91) 



The final reordering is possible because f and A are diagonal and therefore TtFr is diagonal. This equation 
is easily solved for the elements fff to give 



(W)«=I§*=X*/(rtFT)« 



(92) 



The other elements of L are gotten from these values and the ratios g ik already determined. 

The actual motions in the cartesian system can now be obtained by combining eq (39) with the inverse of 
eq (60) or eq (77). For the stretching motions we have 



622 



X = M-'I^M" 



L,Q r 





(93) 



where L r arises from the secular equation containing H r . Q r is the normal coordinate vector and the 
x-component of the center of mass has been set to zero. A similar equation results for the bending modes. 
This is 



Y=A" , ¥=A- 








(94) 



where the z-component of the angular momentum and the y-component of the center of mass have been set 
to zero. 
This completes the frequency analysis. In the next section we will consider the partition functions. 

2.4. Partition Functions. 



Herschbach et. al. [13] have shown how to express the classical partition function for polyatomic 
molecules in terms of local properties. We shall use their method because it allows for cancellations of con- 
siderable portions of the partition functions of the complex and reactants when their ratios are evaluated in 
the rate constant expression, eq (3). We begin the discussion with the classical partition function for a linear 
polyatomic molecule. This is (see Herzberg [14], pp. 502-509), 



q cl = geV(2irMkT/h 2 ) 3/2 (kTI(ahcB)) n (kTlfahc)) 



(95) 



where q e is the electronic partition function, V is the volume, M is the total mass of the molecule, c is the 
velocity of light, o>, is the frequency of the Tth vibrational mode in cm" 1 (a\ = vjc), N is the number of atoms 
in the molecule, B is the rotational constant; B = hl(8-n 2 cl), where /is the moment of inertia of the molecule; 
a is the symmetry number which is the number of indistinguishable positions into which the molecule can be 
turned by simple rigid rotations. For linear molecules a= 1 or 2. Equation (95) neglects nuclear spins, anhar- 
monicity, and non-rigidity of the molecule. Let us rewrite eq (95) in terms of / and Ui — hvJVT. It becomes 



q e ,=q.V4iro- 1 (2T r kTh- 2 ) i/2 M 3/ 



3JV-5 

i n u. 



(96) 



q c i can also be written in the form 



qd—q,o~ x Z II A; 



(97) 



where 



Z=\ \e- vnT dx 1 dz n 

A a = h(2TrM a kT)-» 



(98a) 
(98b) 



Z is the so-called configuration integral, V is the potential energy, and x v y v z x , x N , >V, z N are the 

cartesian coordinates of each of the N atoms. Eliminating q c , between eqs (96) and (97) gives 



z = v^ivkTh-Y^-^M 3 ' 1 1 n M- a 3n 3 ri 5 uj 1 



(99) 



Consider now the matrix H=GF defined by eq (84). A theorem of matrix algebra states that the determi- 
nant of H equals the product of its eigenvalues (see Hohn [15], p. 283). There is also a theorem (Hohn, p. 65), 
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stating that the determinant of product of two matrices equals the products of the determinants of the 
matrices in the product. Consequently, 

3JV-S 3/V-S 

|H|*=|G|*|F|*= n X* = (2TrkT/h) 3N -* n u, (100) 

Solving for the product over u t gives 

3JV-5 

II ul l = (2TckT/h) 3N -* | G | -* | F | -w (101) 

Inserting eq (101) into (99) yields 

Z = F4ir(2xJfcr)*< 3 "- 5) M 3 ' 2 1 U. Ml 3 ' 2 | G | "* | F | "* (102) 

This can be rearranged to give 

| F | » (2%kTy w3N ' S) Z = V4tcM*' 2 I n Ml 3n | G | -* = J„ (103) 

The left side of eq (103) does not involve the masses, while the right side does not contain force constants. 
Therefore, the quantity denoted by J N does not depend on either the force constants or the masses, but must 
depend only on geometrical parameters. Herschbach et. al. [13] have shown that for linear molecules 

J H = F4t n 1 R% u (104) 

where i? m>/ is the equilibrium distance between mass M t and M,*,. For a general linear molecule, the 
classical partition function per unit volume can now be written 

Qc = qJV = V 1 q.<j- l M2-KkT)^ 3N -* |F|-* UK 3 (105) 

This form of the partition function is suitable for the reactant molecules. 

Let us now consider the partition function for the complex. Using eq (96), we have 

kThr 1 q* = kTh~ l q.o- 1 VA*{2-KkTh-Y n M 312 1 'ft 6 i*; 1 (106) 

Note that the product is over 3iV— 6; i.e., one less vibration than in a stable linear molecule. Consider next 
the quantity 

.3/V-6 v 31V-6 3/V-6 

( 1^ ul l ) kThr 1 = (kT/h) 3N -> n pj 1 = (kT/h) 3 ™ n 2ttX7* 

3W-6 

= (kmy*-* (2ir) 3 "- 6 x*a n x;*x-* 

= (2^kT/h) 3N - i p*|F|-w|G|"W (107) 

where X* is the negative eigenvalue and v* is the associated imaginary frequency; eq (100) has been used. 
Using eq (107) in (106) gives 

kTh- 1 Q* = kTh- 1 q*V- 1 = V- 1 q.a- 1 [V^M 312 IU M; 3 ' 2 \G\-^*\¥\-^2TkT) 3N - i/:i 

-h- 3 " II Ml! 2 

a 

= V 1 q.<T x J*v* I F I -v* (2xkT)^ 3N -^ n A; 3 (108) 
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This equation is very similar to eq (105), the partition function of a stable molecule. Note that | F| _Vi will be 
imaginary for the complex. 

We can now write down the specific partition functions per unit volume for the four reaction cases shown 
in figure lb. 

Case III. 

Species A = M 2 —M 3 

Species B — M 4 

Species C = M 2 . . ,M 3 . . .M t 

Qa = q tA o A 1 4TrRlA2TrkT)^F^(A 2 A 3 r 3 

Qb = q tB A4 3 

QckT/h = q^oc 1 IttRIRIv* ]F r |-*F"4 (2irife7? (AAA,)" 3 

The matrix F P is the 2x2 one given by eq. (5b), and not the 4X4 used in eq (83). 

Case IVa. 

Species A = M l —M 2 —M 3 

Species B = M 4 

Species C = M x — M 2 . . .M 3 . . ,M 4 

Qa = q tA a A l 4irR 2 as Rl(2TrkTf(F al F b r^l(A 1 A 2 A 3 )- 3 

Q B — g. B Ai 3 

QckT/h = q*o? 4TcRlRlR 2 c i>*F^\F r \-»F<\F? 3 (2**77" (AAA3AJ- 3 

Note that I have included Fa, in Qa and Q c even though it is supposed to be infinite; it will cancel out when 

the ratio QJQa is taken. Also note that the bending force constants appear with twice the power of the 

stretching force constants. This is because of the degeneracy. 

Case IVb. 

Species^ = M 2 —M 3 

Species B = M 4 —M s 

Species C — M 2 . . M 3 . . .M 4 -M s 

Qa = q^o-A l 4T/?U2 7 rA7TTO(A 2 A 3 )- 3 

Qb = qeB tf 47ri?i (2T*r)* ra (A 4 A 5 r 3 

QckT/h = g ^^U7ri?ii?fi?i^lF r |^F;^(27r^ 7 ' 2 (A 2 A 3 A 4 A 5 r 3 ^ 

Case V. 

Species A = M 1 -M 2 -M 3 

Species B = M 4 —M s 

Species C = M 1 -M 2 . . M y . .Af 4 -Af 5 

Qa = qtA OA l ^RlMA2*kTftFJ b r*Fl\,{KKKY 3 

Qb = gW47T/?L(2^r)^<MA 4 A 5 )- 3 

Q c kT/h = g.^ 4x/?L/?OT?L v*F-^\Y r \-^^F-^Fl\Fl\(2-wkTf (A 1 A a A 3 A 4 A s )- 3 

We now have everything for eq (3) except the tunneling correction. This will be taken up in the next 
section. 



2.5. Tunneling Correction 

The one-dimensional Eckart potential function was used to approximate the barrier to quantum 
mechanical tunneling from reactants to products. Three parameters are required for its definition; these are 
shown in figure 8. Its functional form is 
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nx) 1-y (l-yf 



(109) 



where 



y = -e 2wx/L 



A =V 1 -V 2 

B = (F^ + v^y 

L =2ii~2IF)»(Vl* + V;*Y\and 



F = 



a* 2 



evaluated at the maximum in the curve. F is a force constant. Using this potential function, Eckart [16] 
solved the wave equation and obtained the transmission coefficient for a particle with mass m approaching 
the barrier from the left with an energy E. His result is 



KiE,V v V 2 ,F) = 1 - 



coshp^o:! - oc 2 )] + A 
cosh^^a! H-aJ] + A 



(110) 



where A = cosh[27r6] if 5 is real, and A = cos[2tt 1 5 1 ] if 5 is imaginary. The relationships of <x v <x 2 , and 6 to 
the parameters of figure 8, are 



a, = V2(E/C)» 
« 2 = Vz[(E-A)IC]* 
6 = i/ 2 [(fl-C)/CP 
C = h 2 l(8mL 2 ) 



(HI) 



V 





F 
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t 


s "\ 


v ♦ / \ 



X 



FIGURE 8. Eckart potential function. 



Given the transmission coefficient, Johnston [2], pp. 42 and 43, has derived the correction factor T* which is 
the ratio of the quantum barrier crossing rate to the classical crossing rate. His result is 



T* = eV* r j; K(E)e- E "> T dE/kT 



(112) 



where E„ = when V x < V 2 and E „ = V, - V 2 when V x > V 2 
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Let us rewrite this in a more symmetrica] form. We define a new variable € = (E — VJIkT, E — kTe + V v 
Equation (112) becomes 

r* = j"K(e)e-<fc (113) 

where e a = -VJkT when V x < V 2 and e D = -VJkT when V x > V 2 . 

With this substitution, the parameters a l and a 2 become 

a,- = V 2 {k Te/C + Vt/Q* i = 1 & 2 (114) 

From eq (110) we have 

K = K(a 15 a 2 ,5) = KikTe/C,VJC,VJC,B/C) 

But 

B/C = VJC + 2[(VJC)(VJC)]» + VJC (115) 

is a function of VJC and VJC. Therefore 

K = K(e,p,/> ljP2 ) (116) 

where 

p = kT/C, pj = VJC and p 2 = F" 2 /C. 

T* thus depends on three parameters. Furthermore, it is invariant when p t and p 2 are interchanged; i.e., 
T*{p>Pi,pJ — T*(p,p 2 ,p l ). To see this let p[ = p 2 and p 2 = p v From eq (115) we see that 

(B/C)' = ^ + 2(p' lP ^ + P ' 2 =p 2 + 2(p 2Pl )* + Pl = B/C 

Thus, 8' = 8. From eq (114) we have 

oc[ - Vz(pe + p'J* = Wpe + p^ = a 2 

Cl' 2 — «i 

Using these results in eq (110), we get 

Kie,p, P l,p 2 ) = K(e,p,p v pJ = Kie,p,p v p^ (117) 

Suppose that p[ > p 2 ; i.e., V[ > V 2 . Using eq (117), eq (113) becomes 

r*(p» j p> 2 ') = ir p> K(e,p,p^p 2 )e- t de= \~ plp K(e t p,p l ,pje"'dc=r*ip,p l ,p3 

The way Eq (113) was integrated to get T* will be considered later when the computer program is 
discussed. 

In applying this correction, it is assumed that the x coordinate of Eckart's potential lies in the q direction 
discussed earlier. This is that direction at the saddle point in which the potential energy decreases most 
rapidly. It is also the direction of the path of constant total bond order. We therefore use the force constant 
F„ given by eq (22) for the second derivative of the Eckart potential at its maximum. The effective mass for 
tunneling, M„ is the proportionality factor between the kinetic energy and V2Q 3 . We can calculate M, in the 
following way: As far as tunneling is concerned, in the 4 and 5 mass point cases there are effectively 3 
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masses, since the end bonds are supposed to be rigid. Thus, there are only the two variables, R b and R« 
involved. (Bending modes are not considered.) The kinetic energy T T for changes in these two bonds is given 
by 



T, = VflRfGrK 



(118) 



where R is the 2-dimensional vector defined by eq (6), and G r is the 2 X 2 matrix given by eq (86). The 
inverse of this matrix is easily calculated and found to be 



G; 1 = 



Af 2 (M 3 +A£) 



M^Ml 



M' 2 M; (m 2 '+aqm: 



M" 



(119) 



where M' 2 = M Y + M 2 and M[ = M t + M s in the 5 point case. The transformation between R b , R c and q, a 
is given by the matrix U whose value, determined by the BEBO calculation, is given by eq (18). U can be 
used to express T r in terms of q and a. Thus 

T r = 1/2 RtG^R = i/ 2 PfUfG^UP 

The desired quantity M, is simply the matrix element (IW^Uhj. This is 



M, = 



M'AM S + M 4 ' )m 2 -2M 2 'M 4 W +(M 2 ' + AQM 4 V 
(n* + m z )M 



(120) 



where n and m are the bond orders from the BEBO calculation, and M is the total mass of the molecule. 

The bases from which the tunneling parameters V x and V 2 are measured are taken to be the zero point 
energies of the reactants and products, respectively, and not the potential minimums as might be expected. 
The maximum of the potential, on the other hand, is placed at the potential minimum of the complex; i.e., at 
the saddle point. Johnston [2], pp. 190-196, gives reasons for this particular method of using the Eckart 
function for tunneling corrections. 

We finally have everything needed for eq (3). In the next section explicit rate constant expressions will be 
given for the four reaction cases of figure 2. 



2.6. Rote Constant Expressions 

The rate constant expression eq (3) is not quite complete. It should be multiplied by the number of 
equivalent H atoms on the molecule being attacked. Let us call this factor the chemical multiplicity, <?<*. For 
example, there are 6 identical reaction paths for H abstraction of the 6 terminal H atoms on propane, and 2 
paths for abstraction of the 2 central H atoms. Thus a eh = 6 in the first case, and 2 in the second. With this 
factor added, the rate constants for the four cases shown in figure 2 are 

Case III. M a -M 3 + M t -M s . . .M 3 . , .M A 

k = sn\ (2xir) 3 ' 2 (ry 

Case IVa. M t -M 2 -M 3 +M A -M 1 -M s . . .M a . . ,M 4 

k = S/vJfV^" 1 VirkT) 3 " (r.TD 2 (r 2 T 2 (121) 

Case IVb. M S -M 3 + M 4 -M S -M P . .M 3 . - .M.-M, 

k = SiF^Ftj- 1 (2irkTy» (Tj^r)" 
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Case V. M 1 -M 2 -M 3 +M 4 -Af 5 -M 1 -M 2 . . M y . .M 4 -M 5 

The common factor in all these expressions is 

S = a „ qtC 0a<Jb R2bR ' v* F ^' r * T* e' v ' ,kT 
q^q, B o c Rl, |F r |*r? 

The calculated factors in S are: 

1) R b and R c ; these are calculated from n and m through Pauling's relation, eq (9). 

2) v* is the imaginary frequency obtained from the vibration analysis for the asymmetric stretch. 

3) |F,] is the determinant of the matrix given by eq (30). It is negative. 

4) r* is the quantum correction factor for the symmetric stretching frequency obtained from the vibra- 
tional analysis. 

5) r* is the tunneling correction factor obtained in section 2.5. 

6) V* is the saddle point potential energy given by the BEBO calculation. 
Other calculated factors are: 

1) iv, is the bending force constant given by eq (33). 

2) F+ 2 and F^ 4 are the bending force constants given by eq (34). 

3) The quantum correction factors T*, T*, T*for the bending modes come from the frequency analysis via 
eq(2). 

This concludes the theoretical part of this discussion. The next section contains a brief discussion of the 
computer program which was written to implement the rate constant calculations. This will be followed by 
instructions on how to use it. 



3. Computer Implementation of BEBO 

The computer program consists of a main section and six subroutines. It is written in an enhanced form of 
BASIC for use on a Hewlett-Packard 9845A computer. 

3.1. Description of the Main Program 

The main program begins by reading the following data: 

1) RunidS 

This is a string variable having up to 79 alphanumeric characters to be used for the run identification. 

2) Opt(M),M=lJ 

These are flags which provide a series of available options. These will be described in detail in the instruc- 
tion section. 

3) Ntemp 

This is the number of temperature values at which the rate constant is to be evaluated. A maximum of 16 
values will be allowed. 

4) Trnin/Tmax 

The minimum and maximum temperature values desired. The reciprocal temperature scale is divided 
into Ntemp — 1 equal intervals and the temperature evaluated from the reciprocal values. This gives a better 
distribution on an Arrhenius plot than if the temperature scale were divided into equal intervals. 

5) M1,M2,M3,M4,M5 

These are the five mass point values determined according to the rules given in section 2.1. 

6) Ras,Rbs,Rcs,Rds 

These are the equilibrium bond distances for single bonds. 

7) Ebs,Ecs,P,Q 
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The first two parameters are the electronic energies for single bonds b and c\ the last two are the BEBO 
parameters obtained from eq(12). 

8) Rts,Ets,Beta 

These are the bond distance, bond energy, and Morse parameter ft for the triplet interaction. 

9) Fbs,Fcs,Fpsi2s,Fpsi4s 

These are the stretching force constants for single b and c bonds, and the bending force constants about 
the M2 and M4 masses. 

10) Sa,Sb,Sc 

These are the partition function symmetries o A , o B , and o . 

11) Schem 

This is the chemical multiplicity o c a- 

12) Sea,Seb,Sec 

These are the electronic degeneracies q eA , q, B , and q eC . 

The program next prints out this input data to provide an easily read record and a check of the numbers. 

After these preliminaries, the program then determines the saddle point position. This is done by an 
iterative procedure; n is initially set to 0.5; then the potential energy Pis calculated according to eq (14) 
along with its first and second derivatives, Vn and Van, with respect to n. The subroutine Trpl is used to 
calculate the triplet part of V. A new n is estimated by the Newton, Raphson method from the formula ri = n 
—Vn/Vnn. The process is repeated using n and continued until covergence is obtained. This yields a value 
of n which makes Vn zero; this will correspond to the desired maximum in V. (I have not investigated the 
conditions for which a maximum is expected or if there could be more than one maximum.) 

Having obtained the value of n for the saddle point, the program calculates the stretching force constant 
matrix Fr given by eq (30), its determinant, and the saddle values of Rb and Re from Pauling's relation eq 
(9). It then evaluates the mass to be used for tunneling from a somewhat rearranged eq (120). Next, the 2 X 
2 matrix Gr is calculated from eq (86). This is then combined with Fr to form Hr, and the stretching frequen- 
cies obtained by solving the resulting quadratic secular equation. The bending frequencies are next deter- 
mined through the matrices F(eq (35)) and G (eq (87)). The sizes of these matrices will depend on the type of 
reaction. For the three mass point model there is only one element and thus a linear secular equation with 
one bending frequency. The two four point models require solving a quadratic secular equation for two fre- 
quencies. The five point model uses the subroutine Cubic to solve the cubic secular equation for three fre- 
quencies. The subroutine Normod then calculates the matrix for the normal coordinate transformation of 
the stretching modes. 

At this point, the program prints out a number of properties of the complex. This will be discussed in 
detail in the instruction section. 

The rate constants are then evaluated from eqs (121) at the different temperatures. The activation energy 
is gotten by numerically differentiating the logarithm of the rate constant by means of suitable finite 
difference formulas. Subroutine Fit is a least-squares routine which is used to fit Arrhenius equations 
through the calculated points. The program concludes with subroutine Pltk which draws an Arrhenius plot 
of the results. 



3.2. Discussion of Subroutine Tun 

The only subroutine worth discussing is Tunl, the routine for evaluating the integral of eq (112) for the 
tunneling correction factor T*. Johnston and Heicklen [17] calculated this integral numerically by an 
unspecified method for a range of input parameter values. The three input parameters which they used were 
hv'/kT, where v* - (-F/m)' 4 /(2x), 2-kV x \{Kv% and ItV^Kv*). Their results are in the form of a table. The 
method used in the present program is a modified 6-point Gaussain quadrature formula based on Legendre 
polynomials (see Abramowitz and Stegun [18]). This was used even though the nature of the integral sug- 
gests using a formula based on Laguerre polynomials. Neither of these formulas was satisfactory for the 
whole range of parameter values given by Johnston and Heicklen, so a modification of the first method was 
developed. It was based on the following ideas: When e gets large, the transmission approaches unity. The 

630 



idea is to use the Gaussian formula for that part of the integral where K(e) < 1. After K(e) has gotten suffi- 
ciently close to unity, the remainder of the integral can be evaluated analytically; i.e., if Kfe)* 1 for e>e», 
then 

The problem is to estimate e 6 . Let us examine eqs (1 14) as e — oo . We get a,— l A£*, where £ = kTe/C. From 
eq (110), we have 

K - l -(l + AX**exp(27r!%) + A)" 1 = K fl 

We can set K 6 to some arbitrary value close to unity and solve this equation for e and then £ and then i b 
which will be our cutoff point. The result is 

e b = C i(2ir)~Hn [2(1 + A)/(l - K k )]} ^ftT)" 1 

It turns out that this value is not entirely satisfactory and subtracting from this the average value of V x and 
V 2 works better. Also it can happen that e b as calculated from this formula can be very large when K is close 
to unity. Thus, exp( — e b ) will be very small. There is no point in using a value for e as the upper bound to the 
Gaussian formula if the integrand at this point is negligible because of the exponential factor. Thus e b was 
kept below a certain fixed value e mox . This yielded two parameters, K<, and e max which were adjusted to 
minimize the squares of the differences between the results of this method and the results of Johnston and 
Heicklen. The differences averaged 1.3 percent with only two value differing by as much as 6 percent Such 
accurancy should be quite adequate for the rate constant calculations. 

4. INSTRUCTIONS FOR USING BEBO 

4.1. Input 

It will be assumed that the reader is familiar with the general operation and command system of the 
HP9845A. The program lines 5000 to 5240 contain a series of DATA statements which hold the input data. 
As an example, data for the ethane plus methyl radical reaction is contained in these statements. The 
general nature of the input has been discussed briefly in the last section; here this is considered in more 
detail. 

1) RunidS is a string variable containing identifying information; 79 characters can be used. 

2) Opt(M),M= 1,7 are flags for the following options: 

Opt(0): This picks out the version of the triplet function V t ; these different forms of V, will be discussed in 
the Appendix. 

Opf(l): As mentioned earlier, the activation energy Eact at any temperature is obtained by numerically 
differentiating the logarithm of the rate constants. This is done in either of two ways. The more accurate 
method evaluates the rate constant three times at each temperature; at the particular point and slightly 
above and below the point The derivative is then estimated from a 3 point finite difference formula. This is 
automatically the method used when only a single temperature point is requested. The second, less accurate, 
but faster method uses the rate constants calculated at Ntemp(see last section or below) points and uses a 5 
point difference formula for the derivative. The more points requested and the narrower the temperature 
range, the more accurate is this method. The value of Opt(l) determines which of these methods will be 
used. Thus, when 

Opt(l)= 1, 5 point difference formula used to get Eact (fastest method). 
Opt(l) = 2, 3 point difference formula used to get Eact (most accurate method). 
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Opt(2): When 
Opt(2)= 1, the lnatural logarithm of the rate constant is calculated. 
0pt(2)=2, the logarithm, base 10 of the rate constant is calculated. 

OpiSy. When 
Opt(3)=l, the cathode ray tube is used for the printout In this mode, execution of the program pauses 
before the Arrhenius plot is produced, and before the caption to the plot is generated. In each case execu- 
tion can be resumed by pressing the "cont" key. 
Opt(3)=0, the internal printer is used for the output. 

Opt(4): When 
Opt(4)= 1, the rate constant is in cm 3 /mole-s. 
Opt(4)=2, the rate constant is in cm 3 /molecule-s. 
Opt(4)= 3, the rate constant is in liters/ mole-s. 
Opt(4)=4, the rate constant is in liters/molecule-s. 

Opt(5): Not used. 

Opt(6): When 
Opt(6)=0, the Eckart tunneling correction is not applied. It will automatically not be applied if the zero 
point energy of the reactants is greater than the potential energy V* of the saddle point. 
Opt(6)= 1, the tunneling correction is applied. 

Opt(7): When 
Opt(7)=3, the three parameter Arrhenius type equation, AT" 'c-Exact-mr is fit to the calculated rate constant 
values. 

Opt(7)=2, the standard two parameter Arrhenius equation Ae - £act * RT is fit to the calculated rate constant 
values. 

3) Ntemp is the number of temperature values (up to 16) at which the rate constant is to be evaluated. Use 
the absolute temperature scale. 

4) Tmin, Tmax are the minimum and maximum temperature values to be used. If Ntemp= 1, then only one 
temperature value should be entered on this line. 

5) Afl,Af2,Af3,Af4,A/5 are the five mass point values determined by the rules on page 5. For 3 point models 
set Ml and Af5 to zero. The 4 point models will have either Ml or M5 equal to zero. Atomic mass units are to 
be used. 

6) Ras,Rbs, Res, Rds are the single bond distances in A. For 3 point models set Ras and Rds to zero. For 4 
point models, set either Ras or Rds to zero. 

7) Ebs,Ecs,P,Q; the first two parameters are the electronic energies for single bonds in kcal/mole. The 
quantity normally available is the bond dissociation energy DH° which is defined as the enthalpy change in 
the process in which one mole of the bond of interest is broken, with reactants and products being in their 
standard states as ideal gases at 1 atm and 25 "C. This is not the energy we want. The desired energy E is 
shown in figure 9, which illustrates the energy relationships involved in the removal of an H atom from some 
group A. Z A - H and Z A . are the zero point energies for the reactant and molecular product, and Hj»_ w , H*., 
and H«. are enthalpies of the speices A-H, A-, and H-, respectively. In general, a particular enthalpy is the 
sum of the translational, rotational, vibrational, and PV contributions. We have 

IF = iF( trans) + /F(rot) + /F(vib) + PV 

By examing figure 9 it is easy to derive the relationship between E and DH°. It is 

£ = DIT + (m-H - ft.) + (Z A . H - Z A .) - W H . (122) 

632 



T T 
A- 



- ^1 



77 



DII 



u° 



1 






"A-H 



E 



I n 



2 A-H 

ixx 11 




FIGURE 9. Bond energy relationships. 



The second term is 

IF a -h -HJ. = JE^trans) - HI (trans) + /fl-«(rot) - i£.(rot) + fl^vib) 

/£_„ - m. = i£-«(trans) - /ft. 

Assuming equipartition of energy, the translational and rotational enthalpies will be the same and the dif- 
ference in vibrational enthalpies will normally be negligible. Thus, the second term in eq (122) can be 
neglected. The last term HI. = EJ,. + PV = 3RT/2 + RT, where 3RTI2 is the translational energy of the H 
atom and RT is PV for an ideal gas. Thus, eq (122) becomes 

E= DH° = (Z A .„-Z A . )-5R 772 

As an example, consider the process CH 3 ~H-*CH 3 + H '• . To estimate the difference in zero point energies 
between CH 3 —H and CH 3 -, I have assumed that one C*H stretch of 3100 cm" 1 and two H-C-H bends of 1450 
cm" 1 have been lost in going from A-H to A- and H-.This corresponds to a zero point energy difference of 
8.575 kcal. For cases like this, the bond energy will be 

Ecs = DH + 8.575 -5i?T 2 , 8 /2 = DK + 7.095 kcal 

The zero point energy difference for other types of bonds can probably be satisfactorily estimated in a simi- 
lar manner. HavingobtainedfEbsandEcsin this manner we can calculate P and Q fromeq(12). 

8) Rt$,Et$,Beta are the triplet interaction parameters in A kcal and A" 1 , respectively. I have been using the 
values given in Johnston [1966], table 11-1. 
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9) Fbs,Fcs,Fpsi2s, Fpstts are the single bond force constants. The first two are the stretching constants in 
dynes/cm; the second two are bending force constants in dyne-cm. In the 3 mass point case, both the bend- 
ing force constants are set to zero. For 4 point models, only one of the bending force constants will have a 
value of zero. 

10) Sa,Sb,Sc are the partition function symmetries for A-H, B • , and A • • H • • B, respectively. 

1 1) Schem is the chemical multiplicity. 

12) Sea,Seb,Sec are the electronic degeneracies for A-H, B-, and A- -H- 'B. Sea will normally have the 
value one. Since B ■ and A • • H • • B each have an unpaired electron, Seb and Sec will normally have the value 
two. 

4.2. Output 

BEBO first prints out the input data. It then the following properties of the complex: 

1) The potential energy of activation V* in kcal/mole. 

2) The bond orders n and m of the b and c bonds. 

3) The bond distances Rb and Re in A. 

4) The force constant in the q direction in dynes/cm and the angle q makes with the Rb axis on a contour 
plot like figure 3. 

5) The force constant in the o in dynes/cm, and the angle to the Rb axis. 

6) The force constant in the unstable normal mode direction in dynes/cm, and the angle to the Rb axis. 

7) The force constant in the stable normal mode direction in dynes/ cm, and the angle to the Rb axis. Note 
that the normal mode directions are usually not orthogonal. 

8) The stretching force constant matrix Ft in dynes/cm. 

9) The equations for transforming back and forth between the normal mode and valence bond coordinates. 

10) The bending force constants in dyne-cm. 

11) The two stretching frequencies in cm" 1 . 

12) The one to three bending frequencies in cm" 1 . 

13) The zero point energy of the complex in kcal/mole. 

14) The zero point energy of the reactants in kcal/mole. 

15) The zero point energy of the products in kcal/mole. 

16) The Eckart potential function parameters V\ and V2 in kcal/mole. 

17) The reduced mass for tunneling M,= Mrho. 

18) The second two of Johnston and Heicklen's tunneling parameters (see section 2.5). 

The program then prints out the rate constants as a function of temperature. Also given at each temperature 
is the logarithm of the rate constant, the logarithm of the Arrhenius preexponential factor, the activation 
energy, the tunneling correction factor, and the first of Johnston and Heicklen's tunneling parameters. 
Since the tunneling algorithm has not been checked outside the parameter ranges used by Johnston and 
Heicklen, their parameters values are listed to make sure that they are within the proper ranges. The limits 
areAl and/12 = to20, and U* = to 16. 

Finally, there are listed the differences between the calculated values of the logarithm of the rate constant 
and the values predicted by the least squares fitted Arrhenius equation. This fitted curve is shown by the 
dotted line on the Arrhenius plot. The fitted Arrhenius parameters are given in the caption ot the plot. On 
the next two pages there is a sample output for the ethane and methyl radical reaction. 
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6. APPENDIX: Various Triplet Functions 

The subroutine Trpl is able to provide three different triplet functions which are selected according to the 
value of flag Opt(0). They are as follows: 

Opt(0)=0: This is the modified Sato triplet function with a small portion neglected. Instead of Eq. (13), V, 
= E t5 is used. This simpler formula seems to have been used in the days of mechanical desk calculators. This 
option is useful when attempts are being made to reproduce the results of earlier workers. 

Opt(0)=l:Eq.(13)isused. 

Opt(0)=2: Arthur eL al. [19] have developed a triplet energy formula by fitting a function to the H-H 
triplet potential energy values given by Hirschfelder and Linnett [20]. Their formula is 

V t = 5.873E ts e- 1 - 747 ^ li >>*« c ){li{R b +R <: )y^ 
They claim better results in certain cases when this function is used. 



635 



BEBO- Cal cuUt ions 

P u n 1 d e n t i f i c at i on : 

CH3CH2-H + CH3 = CH3CH2 + H-CH3 

Op t i oni Us ed i n Ca 1 c u 1 at i oni : 

Modified Sato triplet function. 

F i >.} e - p o i n t. d i f f e r e n c e f o r m u 1 as u s e d to get ac 1 i •• a t i o n e n e r g y . 

B a 3 e 10 1 o g a r i t h m o f t h e r ate c o n s t an t . 

Rate constant units in 1 i t er-s-'mol e-s. 

I'lii s e s , i n ft t o m i c M as s U n its. 

Ml= 17.0510, M2- 12.0UQ, M3= 1.00S9, M4= 12.0118, M5 = 3.02 
S i n g 1 e B o n d B i s t an c e s , i n fl n q s t r o m s . 

R as = 1 . 52600, Rbs= 1 . 09000, Re s= 1 . 09000, Rds=1.09008 
S i n g 1 e E o n d E n e r g i e s o f C e ri t e r E o n d s , in k cal; a 1 s o p & q P an am e t e r 

Ebs=105.100, Ecs=l 1 1 . 100, p= 1.085, q= 1.093 
Sinqle Eond Energy, Distance, '& horse Parameter for Triplet Interac 

Ets = S4.4O0, Rts = 1.54000, Bet a= 1.4250 
Single Bo n d S t r etc h i n g F o r c e C o n s t a n t s in d y n e s •-•' c rn f o r C enter B o n d s 

Fbs=4. 79000E+05, Fcs=4. 79OO0E+O5 
Single Bo n d E e n d i n g F o r c e C o n s t a n t s i n d y n e - c m f o r u t e r M a s s e s . 

Fpsi 2s=9. 14337E-12, Fpsi 4s=5. 46530E-12 
Partition Function Symmetry Numbers for Species ft, B, & C. 

Sft=l, SB=1, SC=l' 
Che m i c a 1 flu 1 t i p 1 i c i t y . 

Sc hem=6 
Electronic Degeneracy for Species fl, B, & C. 

3eft=l, 3eB==2, SeC=2 
Pauling's Bond-Order Parameter. 

Lamda=0. 28 00 

Properties of Complex. 

Potential Energy of Activation: V= 14.596 kcal 

Bond Order Parameters: N=0.5804, M=0.4196 

Bond Distances for Center Bonds: Rb=l. 24235, Rc=l. 33313 Angstroms 

Force Const, in Rho Direction." Frhc= -8.25341E+04 dynes -cm: Angle 

Force Const, in Sigma Direction: Fs i qma=2 . 90227E+05 dynes. -'cm: Angle 

Force Const, in Qr n.m. direction: Fqr = -7 . 53301 E t-04 dynes 'cm: Angle 

Force Const, in Qs n.m. direction: Fqs- 2.21465E+05 dynes c r.i : Angle 

F M a t r i \' f o r S t r e t c he = i n djni e s •■-' c m 
1.62224E+05 1.77020E+05 
1 . 77020E+05 4. 54138E+04 

Norma ! Coord i n at e Tr ans for mat i ons 

0s= 3.2115 Rb + 3.0558 Re Rb= 0.2649 Qs + -0.9779 Qr 

0r= -0.1525 Rb + 0.8279 Pc Rc= 0.0483 Qs + 1.027 7 Qr 

F H a t »•• i y. E 1 e m e n t s f or B e n d s ; < t h e s e e q u a 1 t h e b e n d i n g f o r c e c o n s t a n 
Fps i 2=5 . 30926 E- 1 2 , Fps t 3 = 7 . 5*6 1 3E- 1 3 , F ps i 4=2 . 2935 1 E- 1 2 dyne -c m 

Stretching Frequencies: 1610. 40i, 522.47 m .»•..■«■ nos. 

Bending Frequencies: 1132.37, 148.99, 440.72 waue nos. 

Zero P o int E n e r g y o f C o rn p 1 e •■< = 5.6 6 9 k c a 1 



4 



54. 


, 13 


de 


q 


o=; 


, £7 


de 


q 


■4 6. 


,42 


de 


q 


1 Q . 


.44 


de 


q 



Zero Point Energy of Re act ant s= 7.906 kcal 

Zero Foint Energy of Product ~= 7.478 kcal 

Energy Bai-es for Eel art Tunnel tng: Vl= 6.690, V2=13.117 kcal 

Reduced Mass for Tunneling: Nrho = .9314 

Johns t on ■■' i tunneling parameters: Al= 8.812, A2=17.2^9 
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****** 


**** ** 


*********** 


******** 


* * * * * * * * 


^-■^■^^•^■■^"^■■jt*^ 


******** 


************** ********* 


Temp . , 


Rate 


Const . , R f 


actor, fl 


ct. Ener 


qy, Eckar 


t U* Fac 


t or, D e ' 


-■i at i" on from Fi t 


10Q0'-T 


T 


K 


Log(K) 


Log<ft> 


Eac t 


Gam 


u* 


Fit 


6.50 


2000 


1 . 796E+09 


9.254 


1 1 . 724 


22. 5 '? 7 


1.041 


1.200 


-0. 018 


0.73 


1286 


1.030E+03 


•9.013 


11. 179 


18. 626 


1.091 


1 . 367 


0.019 


1 . 06 


947 


S. 1O6E+06 


6. 959 


10.718 


16.294 


1 . 166 


2.534 


0.017 


1 . 33 


750 


1 . 036E+06 


6.015 


10.347 


14.366 


1. 273 


3. 200 


0. 006 


1.61 


621 


1 .396E+05 


5. 145 


10.622 


13.850 


1 . 422 


3 . 867 


-0.006 


1 . 89 


529 


2. 129E+04 


4.328 


9. 727 


13.079 


1.627 


4. 534 


-0.013 


2.17 


462 


3. 577E+03 


3.553 


9.454 


12.462 


1.909 


5.201 


-0.016 


2.44 


409 


6. 501E+02 


2.813 


9. 193 


11.943 


2.300 


5. 867 


-0.011 


2. 72 


367 


1.264E+02 


2. 162 


8. 944 


1 1 . 502 


2.851 


to . 534 


0. 091 


3.90 


•:• n -~, 


2.599E+01 


1.415 


8.716 


11. 135 


3.627 


7.201 


0.020 


**?,*** 


* * * * ************* 


******** 


******** 


********* 


******** 


*********************** 
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Figure 1. Plot, of logar 
the % e m p e r a t u r e . 
ft least squares fit of th 
k =fi* < T-'^n > *EXP < -Ear-r/RT ) , 

Log<fi)=-2.362E+00, n= 3.321, 



t h m o f r at e c o n s t an t 



function of the reciprocal of 
t he flrrheni u -i e x p r e s 5 ion, 



calculated rate constant to 

% he fo 1 1 ou i ng Ma 1 ues for t he par ainst er 

Earr= 8. ?74 



'ields 
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46 



Tpfi 
DIM 



This is 
con stan 
bon d-en 
It i s w 
on a He 
OVERLAP 
READ H, 
DATA 6 . 
f ac=2 
Run 
DIM Opt 
DIM Gam 
COM Del 
READ Ru 
READ Op 
READ Nt 
IF Ntem 
PRINT " 
BEEP 
STOP 
L6 2; IF 
REDIM T 
REDIM E 
IF Ntem 
IF Ntem 
IF Ntem 
IF Ntem 
Rtmi n=1 
Rtmax=1 
IF Ntem 
FOR L=1 
Temp ( L] 
Tpkt kca 
Ktkca I ( 
NEXT L 
De L rt = D 
READ M1 
IF M1 >0 
IF M5>0 
READ Re 
IF Ras> 
IF Rds> 
READ Eb 
READ Rt 
READ Fb 

READ Sa 
READ Sc 
READ Se 
PRINTER 



o Listing of the computer pro 
ts of hydrogen atom transfer r 
e rgy-bond-orde r method (custom 
ritten in an enhanced form of 
wlett-Packerd 9845A computer. 



gram for calculating the rate 
eactions according to the 
arily referred to as BEBOJ. 
the BASIC Language for use 



Null 



16) 



K.Lam, Pmass , Econv f Pred ,Imax ,C , .iu u i. _ -,„_..„ „ 

6 23 4E-27.1 J 8033 E-1 6 , . 28 , 1 .6 59^2E-24 , 6 . 9461 2E-1 4 f 1 E-9 » 1 5 , 2 . 99776 E1 , 

*PI*K*1 .439651 E13 

1 d $ C 7 9 ] 

[7], Tern p[1 :16) . Tpkt kcaL [1:16], 

prM :16] ,LkqprM :16) f LaprM :16 

ta t Lam,Nbebo,Rbs,Rcs.Rts,Betab 

nid$ [Run identifica 

t[D] ,0pt(1 3 ,0pt(2) ,0pt [3) ,Opt( 

emp INumber of tempi 

p<=16 THEN L62 

E r r o r-03 " 



KtkcaLM ;16) ,Kqpr (1 ;16] 
) ,Ust r (1 :16 ) .Error (1 :16 
.Betac f Beta,Ebs,Ecs f Ets 
t i on . 

4] ,0pt (5) ,Opt [6) ,Opt (7) 
ie ratu re values. 



Eactpr (1 
, Pa rt 1 :3 



Ntemp<5 THEN Opt (1 J=2 
emp(1 :Ntemp) ,Tpktkcal(1 :Nte 
actpr (1 : NtempKGemprM ;Ntam 
p=1 THEN READ Tmax !A 
p=1 THEN Tmin=1 
p=1 THEN Opt [1 ]=2 
p>1 THEN READ Tmin,Tmax IT 
000/Tmax 
000/Tmi n 
p>1 THEN DeLrt=(Rtmax-Rtmin 

TO N t emp 
=1 000/ ( Rtmi n+CL-1 )*Delrt) 
L(L)=Tpfac*Temp[L) 
LJ=Tpktkcal(L]/(2*PI) 

eLrt/1000 
.M2,M3,M4,M5 

THEN Mu1=1/M1 

THEN Mu5=1/M5 
s , Rbs , Res . Rds 
THEN Rhoa=1/Ras 
THEN Rhod=1/Rds 
s , Ecs , P , Q 
s , Ets , Be ta 
s,Fcs,Fpsi2s,Fpsi4s 

, Sb , Sc 
hem 

a . Seb i Sec 
IS 16 



mp ) 

P> r 
si 



he 

]/[ 



, Ktkca L [1 :Ntemp) f Kqp r [1 :Ntemp) 

Lkqpr ( 1 :Ntemp) , Lapr 11 :Ntemp) ,Error(1 :Ntemp 

ngle temperature entered. 

minimum and maximum temperatures. 

Ntemp-1 1 



SAtomic mass units. 



! Angst roms 



lEnergies in kcal/mole. 

IStretching constants in dynes/cm 
! Bending constants in dynes-cm. 
[Partition function symmetries. 
[Chemical multiplicity. 
[Electron ic degeneracy. 






470 
4BD 
490 
500 
510 
520 
530 
540 
550 
560 
570 
580 

590 

BOO 
610 
620 
630 
640 
650 
660 
670 
6B0 
690 
700 

710 
720 
730 
740 
750 
760 
770 
7B0 
790 
BOO 
B10 
820 
830 
840 
B50 
860 
B70 
B80 
890 
900 
910 
920 
930 



L7 2: 

PRI 
PRI 
PRI 
PRI 
PRI 
PRI 
IF 
IF 
IF 
IF 
L6 4: 

IF y 

IF V 
| p 

IF 

IF 

IF 

IF 

IF 

PRI 

PRI 

PRI 

Fo rm 

D 

PRI 

PRI 

Fo rm 
PRI 
PRI 

Fo rm 
PRI 
PRI 

Fo rm 
PRI 
PRI 

Fo rm 
PRI 
PRI 

Form 
PRI 
PRI 

Fo rm 
PRI 
PRI 

Fo rm 
PRI 
PRI 



NT 
NT 
NT 
NT 
NT 
NT 



IF 0pt(3]=0 THEN PRINTER IS 

u.^^.^^^^^^^^^,^ BEBO Calculations " 

T Hont i f i r> a -t- 4 n n • H 



"R 
Ru 



it * 
"D 



un Identification:" 
ni d$ 



Opt [0 
OptfO 
Opt(0 
Opt (0 
IF 



]=1 
]=2 
]=3 
Opt 



ons U 
THEN 
THEN 
THEN 
THEN 

(1 3=1 



sed i n 
PRINT 
PRINT 
PRINT 
PRINT 
THEN P 



* ## * 

Ca Ic 
"Sim 
"Mod 
"Two 
"Thr 
RINT 



illations;" 

plified, modified Sato tripLet function." 
ified Sato triplet function." 

-parameter Arthur et.al., triplet function." 
ee-parameter Arthur et.al., style triplet function." 
Five-point difference formulas used to get activation 



0p t t(1)=2 THEN PRINT "Three-point difference formuLas used to get the activation 



ene rg 
ene rg 



Opt [2 
Opt (2 
Opt (4 
Opt 4 
OptU 
Opt (4 
Opt [6 
NT "* 
NT "M 
NT US 
at31 : 

NT "S 
NT US 
at32: 
NT »S 
NT US 
at33 : 
NT "S 
NT US 
at34: 
NT "S 
NT US 
a t3 5 : 
NT "S 
NT US 
at36 : 
NT "P 
NT US 
at3 7 : 
NT "C 
NT US 
at3B : 
NT »E 
NT US 



)=1 

l=? 

3=2 
]=3 

I -8 

### 

ass 

ING 

I 

i no 
ING 

I 
i ng 
ING 

I 
i ng 
ING 

I 
i ng 
ING 

I 

i "9 

ING 

I 

a rt 

ING 

I 

hem 

ING 

I 

I ec 

ING 



THEN 
THEN 
THEN 
THEN 
THEN 
THEN 
THEN 

aje $ a$t sfc <$e 

es , i 
Fo rm 
MAGE 

Le Bo 

Fo rm 
MAGE 
L a Bo 

Fo rm 
MAGE 
L e Bo 

Fo rm 
MAGE 
I e Bo 

Fo rm 
MAGE 
I e Bo 

Fo rm 
MAGE 
i t i o n 

Fo rm 
MAGE 
i ca L 

Fo rm 
MAGE 
troni 

Fo rm 



PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 

$ ape jJk ^t 3$t a$t ajg 

n A t o m i 
at31 :M1 
» M1=" 



"Nat 

"Bas 

"Rat 

"Rat 

"Rat 

"Rat 

"Tun 
# **# 

c Ma 
,M2, 
,3D. 



ural Logarithm of the rate constant." 
e 10 logarithm of the rate constant." 

constant units in cc/mole-s." 

constant units in cc/mo L ecu Le-s . " 

constant units in I i te rs/mo I e-s . " 

constant units in I i ters/mo lecu I e-s . " 
noting correction not applied." 

ft**********************************!! 



e 
e 

e 
e 

s s Units 



M3,M4»M5 

4D,», M2=",3D.4D, », 



M3=»,3D.4D, », M4='»,3D.4D,", M5=",3D.4 



nd Distances, in Angstroms." 
at32;Ras ,Rbs , Res ,Rds 
" Ras=» r Z.5D ?> " J[> ' Rbs=",Z.5D,», 

Bonds , in 



Rcs=»,Z.5D, ", 
kca L ; a lso p & 



Rds=",Z.5D 
q Parameters." 



nd Energies of Center 

at3 3jEbs,Ecs,P,Q 

" Ebs=»,3D.3D,\ Ecs='',3D.3D » p=",Z.3D,", q=",Z.3D 

nd Energy, Distance, & Morse Parameter for TripLet Interaction. 

at34 ; Ets , R ts, Beta 

" Ets=»,3D.3D " Rts=",Z.5D,», Beta=",Z.4D 

nd Stretching Force Constants in dynes/cm for Center Bonds." 

at35 ; Fbs , Fes 

" Fbs=" ,Z.5DE, ", Fcs=",Z.5DE 

nd Bending Force Constants in dyne-cm for 

at36:Fpsi2s,Fpsi4s 

» Fpsi2s=" Z.5DE " Fpsi4s=",Z.5DE 

Function Symmetry Numbers for Species A, 
at37 :Sa,Sb.Sc K ' 

" sA=» r Z l ft , SB=",Z,", SC=»,Z 
Multiplicity." ' 

at38jSchem 
Schem=",Z 



Outer Masses." 
B, & C." 



c Degeneracy for Species A. B, & C." 
at39;Sea,Seb,Sec 



ON 

© 



940 I 

950 

960 

970 I 

980 

990 

1000 

1010 

1020 

1030 

1040 

1050 

1060 

1070 

1090 

1100 

1101 

1110 

1120 

1130 

1140 

1150 

1160 

1170 

11 BO 

1190 

1200 

1210 

1220 

1230 

1240 

1250 

1260 

1270 

1280 

1290 

1300 

1310 

1320 

1330 

1340 

1350 

1360 

1370 

1380 

1390 

1400 

1410 

1420 

143 



: rmat39 ; IMAGE '» SeA=«»,Z » SeB = » ,Z h » , SeC=»,Z 

PRINT "Pauling's Bond-Order Parameter. 



Svm"schQm*Sec*Sa*Sb/CSBajSeb*Sc) 
Fbs=1 E-16*Fbs/Econv "" 



Force constants are converted to kcal energy units. 
Fcs = 1E-1B*Fcs'/Econv 
Fosi 2s=Fpsi 2s/Econv ^a.-* 

I The saddle point position is calcuLated. 

N=.5 

IcountpO N=.5*NoLd 

IF N>1 THEN N=.5*(Nold+1] 

gXtrTrpLrRbB.Rcs.Rta.EtB.Beta.Lam^.Vt.Vtn.Vtnn.Vtr.Vtrr.OpttD]] 

v=Phq*f 1-N~P ]-Ecs*M*Q+Vt , , 

?p N ABsT[N-NSLd]/NoLd]<Preo1 THEN L1 

Icount=Icount+1 

IF Icount>Imax THEN LI 2 

GOTO L11 _ n . „ 

L12: PRINT "Error-01" 

STOP a.*** 

I Next* the stretching part of the force constant matrix ,s calculated. 

Nsq=N*2 

Msq = M*"2 t 

N2m2=Nsq+Msq 

Nm=N*M , 

Frhn = Vnn*Nm"2/(N2ni2*Lam*2 J , ,..„ n 

FSlgraa2"FbS*N-a+Fcs*M-3 + VtPPj/N2ni2 

DIM FrM:2,1:2) ^ . .., . ,.,„,- 

FrM J ]=[Frho*Msg+Fsigma*Nsq]/N2m2 
E ?l=Fr[2.1 ]=f-Frho+Fsigma J*Nm/N2m2 

Fr 2;2l = [Frho*Nsg + Fsign,a*Msq)/N2m2 

Rb=Rbs-Lam*LOG[NJ 
Rc=Rcs-Lam*LOG[MJ 

Df?=SGN(Dfr)*SQR(ABS[Dfr]] 

Jfac=[Rb*Rc/Rbs)*2 

Rhob=1/Rb 

Rhoc=1/Rc 

Cc=-Nbebo/(1-Nbebo) 



144D Ma = M1+M2 

145D Mx=M3 

1460 Mb=M4+M5 

1461 Maxb=Ma+Mx+Mb 

14^1° S£fi°A4r?T^:r:S? C iLSri M ?3 M 1 X Tg?^ 2+Ma * Mx 3/tMaxb*C1 + Cc^]) 

]}?2 Azr [M ="[Mb + Mx /Maxb' ' ' 
KJ ^r 1,2 =Azr 2 2 =-Mb/Maxb 

1480 t Z M r Jh^1 i = tMa + Mx]/MaXb 

i5oS i N**r%*hr?r**r********r*r***^^ 

JcJn Next, the Gr matrix is calculated TMcs «,<ii u n \ a *•„„ o /i pc *. »-r***,, 

1510 with rigid bonds on the ends. wlU h ° Ld for 3 ' 4 ' &5 atom modeLs 

1530 DIM Grl^^l-iP 6 ° f m ° dBL be1ng USBd is determined. [Note, aLL are Linear) 

1540 ModeL=5 

^ §50 IF [M1=01 OR (M5 = D] THEN Model=4 

1560 IF (M1=0 AND (M5=0) THEN Model=3 

1570 ModeL2=Model-2 aBL J 

1580 Mu2=1/M2 

1590 Mu3=1/M3 

1600 Mu4=1/M4 

on lion g r !Ml=M u ?/M+M1/M2]+Mu3 

°^ 3i§H Gr 1 ,2 =Gr(2,1 ]=-Mu3 

h- 1630 Gr(2,2)=Mu3 + Mu4/H+M5/M4) 

1650 DIM t Hrll Cb 2 n ?-2) eqUenCieS evaLuatB d. 

1660 MAT Hr=Gr*Fr" 

1670 Bh = HrM ,1 )+Hrt2,2) 

1680 Ch = Hr 1 1 ]*Hr 2 2 -HrM 21*HrfP 11 

1690 Dh^SQRlBhV^Ch) Mrl1 ' dJ »rl2 f 1] 

3?PS DIMEv 5 ^!? 1 »Frqs[1 :2] ,Lst1 :2] 

HIS P vs 1 =.5*(Bh+DH 

1720 Evs(2J=.5*(Bh-Dh 

1730 FOR 1=1 TO 2 

1750 ^3f4 1=SGNlEvs[I ^* 6B2 - 4 27*SQR(ABS[Evs[I)))/(2*PI) n n 1/ cm . 

1770 1 nistnrjr**********************************************^ 

i4hR c Be S d 1 PP*I re .H uer )9 1 BS nQ w calculated. 
17B0 Fpsi3=-RB*Rc*Vtr/(Rb+Rc ] 

1830 ! ••........,.,, 

1£2R ,' I hrB .? atam model." ' 

1860 Ff^^)iFoli3 MU2 * Rh0b * 2 + Mu4 * RhOC " 2 + Mu3J!t(Rhob + Rhac ^ 2 
1870 MAT'h=G*F 



*n3 



18BD 
1B9Q 
1900 
1910 
1 920 
1930 
1940 
1950 
1 960 
1 970 
1980 
1990 
2000 
2010 
2020 
2030 
2040 
2050 
2060 
2070 
2080 
2090 
21 00 
2110 
21 20 
21 3 
2140 
21 50 
2160 
2170 
21 80 
21 90 
2200 
2210 
2220 
2230 
2240 
2250 
2260 
2270 
2280 
2290 
2300 
2310 
2320 
2330 
2340 
2350 
2360 
2370 



EvbM )=HM ,1 ) 
GOTO L20 

! 

! Four a torn models. 

La4: IF M1 =0 THEN La4sub 

GM ,1 =My1 *Rhoa-2+Mu3*Rhob-2+Mu2*[Rhoa+Rhob]-2 

G 1 ,2 =p(2,1 ]=-Rhob*( (Rhoa+Rhob)*Mu2+[Rhob+Rhoc)*Mu3] 

G 2,2 =Mu2*Rhob~2+Mu4*Rhoc~2+Mu3*(Rhob+Rhoc] "2 

F 1 ,1 }=Fpsi2=Fpsi2s*N 

F(1 ,2 =Ft2.1 )=& 

F(2, 2 J=Fpsi3 

MAT H=G*F 

GOTO La4end 

La4sub: G 1 1 .1 ) =Mu2*Rhob ~2+Mu4*Rhoc~2+Mu3 * ( R hob + Rhoc ) *2 

G(1 ,2 =G(2,1 )=-Rhoc*[lRhob+Rhoc]*Mu3+(Rhoc+Rhod)*Mu4] 

2 =Mu3*Rhoc-2+Mu5*Rhod~2+Mu4*(Rhoc+Rhod)-2 

1 =Fpsi3 

2)=F(2 f 1 )=0 

2 ) = Fpsi 4 = Fpsi 4s*M 



id J = h p i 
H = G*F 



G[2 

FM 

F(1 

F(2 

MAT 

La4end: Bh = HM ,1 )+H(2.2) 

Ch = HM ,1 ]*H(2,£]-HM ,£]*H(2,1 ) 

Dh=SQR(Bh-2-4 i Ch) 

EvbM ]=.5*(Bh+Dh] 

Evb(2)=.5*(Bh-Dh) 

GOTO L20 



! Five 
La5: G 
G(2,2) 
G 3 3) 
G(1 ,2) 
G(2,3 
G(1 3) 
MAT F = 
F(1 ,1 ] 
F ( 2 2 ) 
F(3,3) 
MAT H = 
DIM Rz 
Rz(0)= 
Rz(0)= 
Rz 1 = 
Rzh) = 
Rz(2)= 
CALL C 



atom mo 
[1 ,1 ]=Mu 
=Mu2*Rho 
=Mu3*Rho 
=G(2,1 )= 
=G(3,2)= 
=G(3 1 )= 
ZER 

= F p s i 2 = F 
= F p s i 3 
- F p s i 4 = F 
G*F 
(0:2) 
H(1 .1 )*H 
RzlO)-H( 
H(3 .1 ]*H 
Rztl ]-H( 
H(1 ,1 )+H 
ub i c ( Rz ( 



1 *Rhoa-2+Mu3*Rhob-2+Mu2*(Rhoa+Rhob)-2 
b*2+Mu4*Rhoc-2+Mu3*(Rhob+Rhoc)-2 
c-2+Mu5*Rhod-2+Mu4*(Rhoc+Rhod) "2 
-Rhob*( (Rhoa+Rhob)*Mu2+TRhob+Rhoc)*Mu3) 
-Rhoc*( [Rhob+Rhoc]*Mu3+[Rhoc+Rhod]*Mu4 3 
Mu3*Rhob*Rhoc 

psi2s*N 

ps i 4s*M 

(2.2]*H(3.3)+H(2.1 ] *H ( 3 2 ) *H (1 .3 ) +H [ 3 .1 ) *H 1 1 . 2 ) *H ( 2 .3 ) 

2.2)*H[3.1 )*H[1 .3]-H[1 ,'f)*H[3.5)*H[2 f 3]-H[3 f 5j*H(2,'f ]*H[1 

[1.3]+H[3,2]*H(2.3)+H(2.1 *H 4,2] 

1 ,-f J*H[2,§)-Ht1 f -f ]*Ht3 f 3]-Hl2 f 5j*Hl3,3] 

C2,2]+hT3.3) 

*) Evb(*) I 



! Bending frequencies. 
L20; FOR 1=1 TO ModeL2 
F rqb (11=6 82.427 *SQR(Evb [I] ]/{2*PI] 



CO 



238D 


NEXT 


2390 


PRINT 


24DQ 


PRINT 


2410 


Fo rma 


2420 


PRINT 


2430 


Fo rma 


2440 


PRINT 


2450 


Fo rma 


2460 


PRINT 


2470 


Fo rma 




DDZ 


2480 


PRINT 


2490 


Fo rma 




DDZ 


2491 


DIM L 


2492 


CALL 


2493 


MAT L 


2494 


DEF F 


2495 


Z=Lq( 
Zr = AT 


2496 


2500 


PRINT 


2510 


Fo rma 




MDD 


2511 


Z=Lq( 
Zs = AT 


251 2 


2520 


PRINT 


2530 


Fo rma 




MDD 


2531 


MAT A 


2532 


! MAT 


2533 


DIM R 


2534 


Rcomp 


2535 


Rcomp 
MAT Z 


2536 


2537 


! MAT 


2540 


PRINT 


2550 


PRINT 


2560 


PRINT 


2600 


PRINT 


2610 


PRINT 


2620 


Fo rma 




4D, 
PRINT 


2630 


26 40 


Fo rma 




4D, 
PRINT 


2650 


2660 


PRINT 


2670 


Fo rma 


2680 


PRINT 


2690 


Fo rma 



"P rope r 

USING F 

t1 : IMAG 

USING F 

t2: IMAG 

USING F 

t3: IMAG 

USING F 

t4: IMAG 

.DD " de 

USING F 

t5: IMAG 

.DD" de 

qi (1 :2,1 

N o r m o d i H 

q=INV(Lq 

NGg[Z]=F 

2,2]/Lq( 

N(Z] 

USING F 

t4a: IMA 

Z.DD," d 

2.1 )/Lq( 

NlZJ 

USING F 

t4b: IMA 

Z.DD," d 

Lg = Az r*L 

PRINT A 

comp L ex ( 

LexTl ]=R 

Lex(2J=R 

eq = Az r*R 

PRINT Z 

"F Matr 

USING " 

USING " 

"Normal 

USING F 

t45: IMA 
II Qr n 

USING F 
t46: IMA 

II Qpll 

11 F Matr 
USING F 

tB: IMAG 
USING F 

t7: IMAG 



ties 
o rma 
E "P 
o rma 
E "B 
o rma 
E "B 
o rma 
E "F 

9" 

rma 
E "F 

9" 
:2), 

p(* 
i ) 
rM . 

1 ,2) 



mp Lex. " 

aL Energy of Activation: V=",DD.3D," kcaL" 

der Parameters: N=",Z.4D,", M=",Z.4D 

Re 

stances for Center Bonds: Rb=",Z.5D,", Rc=",Z.5D," Angstroms" 

o*Econv*1E16 , ATN [ Cc ) *1 80/PI 

onst. in Rho Direction: Frho= ",MZ.5DE," dynes/cm:"," AngLe=",M 

gma*Econv*1E16 , ATN ( -1 /Cc ] *1 80/PI 

onst. in Sigma Direction: Fsigma=",Z.5DE," dynes/cm:"," AngLe=",M 

,1 :2] 

I f Fr(*) f Lqi (*)) 



o rmat4a ; FN 
GE "Force 
eg" 
1,1] 

o rmat4b ; FN 
GE "Force 
eg" 

1:2) ,Zeq(1 
b 
c 

comp Lex 
eq 

i x for S t r 
2X,MZ.5DE, 
2X,MZ.5DE, 
Coord i nat 
o rma t45 ; Lq 
GE " Qs= 

o rma t46 ; Lq 
GE " Qr= 

ix Element 
ormat6:Fps 
E " Fpsi2 
o rma t7 ; -Fr 
E "Stretch 



Gg(Zr)*Econv*1E16 ,Zr*1 80/PI 

Const, in Qr n.m. direction: Fq r=" , MZ . 5DE , " dynes/cm:"," AngLe=", 

Gg(Zs)*Econv*1E16 ,Zs*1 80/PI 

Const, in Qs n.m. direction: Fqs=" , MZ . 5DE , " dynes/cm:"," AngLe=", 



:3) 



etches in dynes/cm" 

2X";Fr(1 ,1 )*Econv*1E16,Fr(1 ,2)*Econv*1E16 

2X";Frt2 f 1)*Econv*1E16 f Fr(2 f 2)*Econv*1E16 

e Transformations" 

i (1 ,1 ) ,Lqi [1 ,2) ,Lq(1 ,1 ) , Lq M .2) 

" f MZ.4D,' 1 Rb + " f MZ.4D," Rc%" Rb= ",MZ.4D," Qs 



1 (2,1 ) ,Lqi (2,21 ,Lq(2,1 ) , Lq(2,2) 
•\MZ.4D," 1 Rb + M ;MZ.4D," Rc r ',' 



Rc= ",MZ.4D," Qs + 



",MZ. 
",MZ. 



s for Bends; [these equal the bending force constants)" 

i2*Econv,Fpsi3*Econv.hpsi4*Econv 

="Z.5DE " Fpsi3=",Z.5DE,", Fpsi 4= " , Z . 5DE , " dyne-cm" 

qs(2) ,Frqs(l ) 

ing Frequencies: ",4D.DD,"i, ",4D.DD," wave nos." 



£ 



27QD 
2710 
272D 
2730 
27 40 
2750 
2760 
2770 
27B0 
2730 
2BQ0 
2810 
2820 
2830 
2840 
2850 
2860 
2870 
2880 
2890 
2900 
2910 
2920 
2930 
2940 
2950 
2960 
2970 
c 8 BO 
2990 
3000 
3010 
3020 
3030 
3 40 
3050 
3060 
3070 
3080 
3090 
3100 
3110 
3120 
3130 
3140 
3150 
3160 
3170 
31.80 
3190 



IF Mode 

IF Mode 

IF Mode 

Fo rmat8 

Zcmp L x = 

FOR 1=1 

Zcmp L x = 

NEXT I 

Zcmp L x = 

PRINT U 

Format9 

PRINT " 
I #$*## 

! The r 
Uni t=1 
ON Opt[ 
L u n i 1 2 : 
GOTO Lu 
Luni t3 : 
GOTO Lu 
Luni t4: 
Luni t1 s 
DEF FNG 
DIM Log 
Npt = 
IF Opt t 
DeLrt=1 
Npt = 1 
L61 ; Ev 
Frqbs=6 
Zreact= 
Ev cs=Fc 
Frqcs=6 
Zprod=F 
IF Mode 
IF [Mod 
Rhocs=1 
Rhods=1 
Evpsi 4s 
F r q p s i 4 
Zprod=Z 
L46: IF 
Rhoas=1 
Rhobs=1 
Evpsi 2s 
F r q p s i 2 
Zreact= 
IF Opt! 
L40: IF 
Zreact= 



, L = 3 , THIN E5 INT USING Formats ;FrqbM ] ,NuLL .NuL L 
-5 thfM EBIBT KilSI? Formats F^b 1 Frqb 2) NuLL 

. ?Mlrl N «m Ri yT USING Format8;Frqb(1 ) .Frqb 2 'Frgb[3] 

: IMAGE "Bending Frequencies: »,4D.DD,», »,4D.DD » " 

i r q s 1 1 j 
TO ModeL2 



4D . DD, " wave nos 



Zcmp L x+2*Frqb ( I ) 



,5*Zcm 

SING F 

: IMAG 
# # $ # $ # 

ate constants 



p Lx*2.85 

ormat9 ;Z 

E "Zero 
***$*### 



4] GOT 

Uni t = 
n i t1 

Uni t=1000 
ni t1 

Uni t = 

I Con 
[Z]=Z* 
ra te ( - 



Lun i t1 
02E23 



85E-3 

cmp L x 

Point Energy of CompLex=" 3D 3D " kcaL" 

******************************************************************** 

are now evaLuated at the different temperatures. 
, Luni t2 , Luni t3 , Luni t4 



6.02E26 
t i n ue 

EXP[-.5*Z)/(1-EXP[-Z)} 
1:1] 



1 )=1 THEN L61 
•987E-3*.05/V 



bs= 
82. 

s*f 
82. 
rqc 
L=3 
eL = 
/Re 
/Rd 
= Fp 
s=6 
pro 
,[M 
/Ra 



Fbs*tMu2+Mu3] 

427*SQR(Evbs)/(2*PI) 

bs/2 

Mu3+Mu4] 

427*SQR[Evos]/(2*PI] 

s/2 

THEN L40 
4) AND (M1O0J THEN 
s 
s 

Si 4 !^lMH 3 * Rhocs ' ,2+M 4 5 *R h °ds 

82.427*SQR(Evpsi4s)/(2*PI] 

(jjHhFpQD&i 4s 

odeL=4) AND ( M1 =0 ) THEN L40 

s 

s 



L46 



2+Mu4*(Rhocs+Rhods]-2) 



/Rbi 

s = 6 B2 .427*SQR(Evpsi2sl/(2*PI ] 
Zre« ' " 
6]=( 
Op- 

Zrei 



'2J 



iact+Frqpsi2s 
=0 THEN TunneL = 
»t [6]=0 THEN LBS 
iact*2.8585E-3 



3220 Fn^I t1 U S IN ? M f ?E m SI 1 ° ;z "act 

If 18 ESi r Sf t l§iMS M ^?..?l B f? Z 5;i3 t EnerflV ° f ""cfnt.= « f DD.3D f « kcal " 

3250 V?=5-! 1 r^ct MAGE " Zer ° *° int Ener 3' of Products-, 3D. 3D," Kca L " 

3260 Tunnel=1 

3270 IF V1>0 THEN LB4 

3280 TunneL=0 

3300 gppJpSSJJp^SSS^RS b f a aC J?er n Sei^?W. ed DBOausB zer ° P"'"* •"-roy of raactants is" 

l§§8 Mf*T v gi3[;i b ?:?ss« z i^ 1 d V2 

3350 ™EKT^uJi NS M Fo1 m a!?3l^o aSBS '" EOkart T »™^S : VI - , DD.3D, «. V2-, 00.30," k cal" 
ddBD Format13: IMAGE "Reduced Mass far Tunnpiinn. Mnh«-n nn * n 
3370 Johnston's tunnel i nn rnr<?«n e S- i I.. I ? L l ng : M rno-",DD.4D 

3390 A1=V1*EcDnv*2*PI/?H*Nul lMrh ° PmaSsll/t2 PI) 

gaag A2==y2*Ecor. V *2*PI/ H*Nu 

3401 TCnneUo 01 AN ° (A * <20) THEN L3Q 

S 3403 G^S T L85 CkQrt tu ™ eLin 9 factor not emulated because A1 or A2 is greater than 20." 

ifi§ S? M ? ya * 5 °^ss N p PRINTER IS 16 

3460 PRINT USING Fo rma t20 : L, NtemD 

3470 Format20; IMAGE "CaLcu atinn ratn mnct- a * h nn m*u ^ 

nAon n Lues." uatcuLating rate const, at » r DD,»"th temperature value out of »,DD •• va 

34B0 Rtemp=l/Temp(U ' 

3490 FOR Nt=-Npt TO Npt 

3500 Rt=Rtemp+Nt*DeLrt 

3510 T=1/Rt 

aiiR JP»<tkcalsub=Tpfac*T 

3 53 Ktkcalsub=Tpktkcalsub/(2*PI] 

3550 ?fSr439)T V/KtkCaLBUbl ^ m * J f ac * C * Fr ^ 2 ) *SQR [ Fbs }/ ( Of r *Uni t ) 

HIS g^S^^ 

3B0D GOTO L30 

2219 !t d4 : IF M1=0 THEN Lb4sub 



§6 20 Kc l = Kf ac*Tpktkca lsub"1 ,5*FDsi2s/[FnRi?*Fnc4Qi 
3630 G4a=FNG(Tf«Frqs(1))MFNG{TpF??b(1??*FNG?ff*Fr 



qb(2] J ) -2 






3B40 

3650 

366D 

3670 

3680 

3690 

3700 

3710 

3720 

3730 

3740 

3750 

3760 

3770 

3781 

3790 

3800 

3810 

3820 

3830 

3840 

3850 

3860 

3870 

3880 

3890 

3900 

3910 

3920 

3930 

3940 

3950 

3960 

3970 

3980 

3990 

4000 

4010 

4020 

4030 

4040 

4050 

4060 

4070 

4080 

4090 

4100 

4110 

4120 

4130 



5 4a nS 4 !^ ( f NG(Tf * Fr q bsl * FNGC Tf*Frqpsi2sr2) 
Kg — b4a *Kc L 
GOTO L30 

G4b = G4b>FN f G(f?*Frqi^ p NG tTf Frqb [1 ) J *FNG?Tf* Frgblal ) ) -2 
Kq=G4b*KcL M 

GOTO L30 

pc 5 A M K?i^^ ac *TP^l59 aLs y bA l- 5 * F Psi2s/(4*PI*Fpsi2*FpBi3*FDsi4] 

|™EN r n??M5ii:iE^4TT f ;Ef?sSi]i:?«ii T, * F? ' b,8,? ^ HB ^ fl * i '' bt3iij 

Kg — G5 *Kc L 

L30; Gam=1 

IF TunneL=0 THEN L87 

l£ 7 l V 5STi5YDSigiS*i^SSsg21C' Frt,0 * Econv * -,E1s ' Mrl,0 * Pn,aSB ' T ' G ' lmJ 

Kq=6.025E23*1E-24*Kq*Gam 

IF Nt=0 THEN Kqtab=Kq 

IF Nt=0 THEN Gamtab=Gam 

Lograte[Nt]=LOG(Kq] 

NEXT Nt 

IF Npt=0 THEN L59 

tp C S^ kl-T 1 -fScM E 7 3 * ( 7, L ?9rS^ [ - 1 3 + LogratB(1])/t2*D B Lpt] 

L59? KQpr?U=Kqtab aP r [ L > =LGT ( Kc l tab *EXP[Eactpp I L]/Ktkca L ( L] I ] 
IF 0ptflj=1 TH|N a Lkqpr(L|=Loqrate(0] 

ystrlL)=Nu*H/(K*Temp(L] ) 

Gampr I U=Gamtab 

NEXT L 

IF 0pt[3)=0 THEN PRINTER IS 

IF Npt=1 THEN L60 

DIM BOM :5].D1 [1 :5] .Di (1 ;5] 

D0MJ=-2. 08^33333333^ 

D0(2)=4 

DO (3] =-3 

D0(4 =1 .33333333333 

DOC 5) =-.25 

D1 M ]=-.25 

D1 [2]=-. 833333333333 

D1 [3)=1 .5 

D1 (4)=-. 5 

D1 (5)=B.33333333333E-2 

Di (1 =8.33333333333E-2 

Di 2 =-.666666666667 

Di { 3 ] =0 

Di (4) = . 666666666667 

Di [5]=-8.33333333333E-2 

MAT Eactpr=ZER 



41 40 FOR J = 1 TO 5 

4150 EactprM }=Eactpr(1 ] +D0 ( J ] *Lkq p r ( J ] 

416 Eactpr t2]=Eactpr(2)+D1 I J] *Lkqpr( J) 

4170 Eactpr lNtemp-1J=EactprtNtemp-1 ]-D1 (J)*Lkqpr ( Ntemp- J+1 ) 

41 BQ Eactpr[Ntemp]=EactprINtemp]-DO(J] *Lkqp r ( Ntemp- J+1) 

4190 FOR Na = 3 TO Ntemp-2 

42 Eactpr ( Na ) =Eact p r ( Na ] +Di ( J ] *Lkq p r ( Na+J-3 ) 

4210 NEXT Na 

4220 NEXT J 

4230 IF 0pt(2]=1 THEN MAT Eact p r = Eact p r * ( -1 .987E-3/De L rt ] 

4240 IF 0pt(2]=2 THEN MAT Eact p r=Eactp r *[ -1 . 987 E-3 *L0G ( 1 ) /De L rt ) 

4250 L60: FOR L=1 TO Ntemp 

4260 Afac = Kgpr(L)*EXP(Eactprl"L]/Ktkcal [Lll 

4270 IF 0ptl2]=1 THEN Lap r I L) =LOG ( Af ac ) 

4280 IF 0pt(2)=2 THEN Lap r [ L] =LGT [ Af ac ) 

4290 NEXT L 

4300 IF (0pt(3]=0) AND [Ntemp>1] THEN PRINT PAGE 

4301 IF Ntemp>1 THEN CALL Fi t ( Opt I 2 ] , Nt amp . Temp [ * J . Lkq pr ( * J , Er ra r ( * J . Pa r ( * 1 . Va r i ance . Opt [ 7 ] ) 
43 10 PRINT ''k*************************************^^ 

4320 PRINT "Temp., Rate Const., A factor, Act. Energy, Eckart U * Factor, Deviation from Fit" 

4330 IF 0pt(2)=1 THEN PRINT "1000/T T K LntK) Ln(A) Eact Gam U* 

Fit" 

4340 IF 0pt(2)=2 THEN PRINT ''1000/T T K Log(K] Log(A) Eact Gam U* 

4350 FOR L=1 TO Ntamp 

4360 PRINT USING Fo rma t1 4 ; 1 000/Temp C L) ,Temp C L ) , Kq p r [ L) f Lkqp r [ L] , Lap r [ L) , Eact p r ( L) , Gamp r ( L] , Us 
t r [ LJ ,Er ror ( L) 

as 4370 NEXT L 

£ 4380 Format14: IMAGE 1 X Z . DD , 2X , 4D, 2X. Z . 3DE, 6 ( 1 X, MDZ .3D ) 

4390 PRINT ''*##******##*#******4**********#***#**###^ 

4400 IF Ntemp=1 THEN L63 

4401 IF 0pt(3)=1 THEN PAUSE 

4410 CALL PLtklOptI*] , Ntemp, Tempt*] ,Lkqpr(*] ,Error(*] ] 

4411 IF 0pt[3)=1 THEN PAUST 



4420 IF 0pt(3]=0 THEN DUMP GRAPHICS 
4430 EXIT GRAPHICS 

4440 PRINT "Figure 1. Plot of Logarithm of rate constant as a function of the reciprocal of" 

4441 PRINT "the temperature." 

4442 PRINT "A Least squares fit of the caLcuLated rate constant to the Arrhenius expression , " 

4443 IF 0f>tl7)=2 THEN PRINT " k=A*EXP ( -Ea r r/RT ) » yieLds the foLLowing values for the parameter 

4444 IF S 6pt(7]=3 THEN PRINT "k=A* [T* n ) *EXP [ -Ea r r/RT ] , yieLds the following values for the par 

ameters:" , 

4445 IF (0pt[7)=2] AND [0pt(2)=1l THEN PRINT USING A r rh en i us1 a ; Pa r [ 1 ) , Pa r [ 3 ] /1 000 

4446 IF (0pt[7]=2) AND [0pt{2)=2] THEN PRINT USING Ar rh en i us2a Pa r ( 1 ) Pa r ( 3 ] /1 000 
4448 Arrheniusla: IMAGE " Ln [ A ] =» , MZ . 3DE , " , Ea r r= " , MDZ . 3D 

4443 Arrhenius2a: IMAGE •' LOG ( A ) = ft , MZ . 3DE '' f Ea r r= f ' , MDZ . 3D 

4453 IF (0pt(7]=3) AND I0pt{2]=1] THEN PRINT USING Ar rh eni us1 b ; Pa r ( 1 ) , Pa r [ 2 ) , Pa r [ 3 ) /1 000 

4454 IF [0pt[7)=3) AND [0pt[2]=2J THEN PRINT USING A r r h en i us2b : Pa r ( 1 ] , Pa r [ 2 ) , Pa r ( 3 ) /1 000 

4457 Arrheniuslb; IMAGE " Ln ( A ] =" , MZ . 3DE , " , n=",MZ.3D,", Ea r r=" , MDZ .3D 

4458 Arrhenius2b: IMAGE » Log ( A ) = " , MZ . 3DE , » , n= " , MZ .3D , " , Ear r=" , MDZ .3D 
4460 L63; PRINTER IS 16 

4470 IF 0pt[5]=0 THEN L82 

4480 ! Position of normal coordinate plot. 

4490 L82: PRINTER IS 16 

4500 SERIAL 

4510 END 



50 DO I Run Identification 

5010 DATA CH3CH2-H + CH3 = CH3CH2 + H-CH3 

5020 I Options: 

5021 ! 0,1 ,2,3,4,5,6 ,7 
5030 DATA 1 ,1 ,2,1 ,3,0 ,1 ,3 

5040 ! Number of temperature values. 

5050 DATA 10 

5060 I Minimum and maximum temperatures, lor a single temperature if Ntemp=1 ) 

5070 DATA 333.3333333,2000 

50B0 ! Masses: M1 . M2 , M3 , M4 , M5 

5090 DATA 17.051,12.011,1.008,12.011,3.024 

5100 I Single bond distances: Ra, Rb, Re, Rd 

<£ 5110 DATA T. 526 ,1 .09,1 .09,1 .09 

S 5120 l Bond energies: Ebs and Ecs: BEB0 parameters: p and q 

5130 DATA 105.1 ,T11 .1 ,1 .0852,1 .0925 

5140 I Triplet single bond distance, bond energy, and Morse parameter: 

5150 I Rts. Ets, and Bete 

5160 DATA 1 .54,04.4,1 .425 

5170 ! Stretching constants Fbs, Fes and bending constants Fpsi2s. Fpsi4s 

5180 DATA 4.79E5,4.79E5,9.14B37E-12,5.4653E-12 

5190 I Partition function symmetries for species A, B, and C 

5200 DATA 1,1,1 

5210 I Chemical multiplicity 

5220 DATA 6 

5230 I Electronic degeneracies for species A, B, and C 

5240 DATA 1 ,2,2 






54QD 
5402 
5404 
5406 
5410 
5412 
5414 
5416 
541 B 
5420 
5422 
5426 
542B 
5430 
5440 
5442 
5444 
5446 
544B 
5450 
5452 
5454 
5456 
545B 
5460 
5462 
5464 
5466 
546B 
5470 
5472 
5474 
5476 
547B 
54B0 
54B2 
54B4 
54B6 
54BB 
5490 
5492 
5494 
5496 
549B 
5500 
5502 



SUB Tr 

! Vers 

! Vers 

I Vers 

DeLrs= 

IF Ver 

M=1-N 

R=Delr 

Rn=Lam 

Rnn=La 

IF (Ve 

P1=5.B 

P2=1 .7 

P3=1 .5 

D=Ets* 

E=EXP[ 

Er=-P2 

En=Er* 

Enn=-P 

E r r=-P 

Rp3=FT 

Vt=D*E 

Vtn=D* 

Vtnn=D 

Vtr=D* 

Vtrr=D 

SUBEXI 

L10 : E 

Er=-Be 

En = Er * 

Enn=-B 

Err=-B 

IF Ver 

Vt=Ets 

Efac=1 

Vtn=Et 

Vtr=Et 

Vtnn=E 

Vtrr=E 

SUBEXI 

L1 5: V 

Vtn=Et 

Vtnn=E 

Vtr=Et 

Vtrr=E 

SUDEND 



?in5n^'§?SA?5 : ^ E 5 S ' Be 5?l. l r ai !l'M' Vt » Vtn . Vt nn f Vtr,Vtrr f V 

nS"? Mo^i 1 - f l e 3'. modl ^ 1ed Sato function ' ' 
ion-1; Modified Sato function. 

Rbs+Rcs-Rts UP et,aL " func tion [2 parameter form). 

sion=2 THEN De L rs=De L rs+Rt r 



e r 8 i o n ] 



s-La 

*M/ 

m*[1 

rsi o 

73 

47 

25 

P1 *B 

-P2* 

*Bet 

Rn 

2*Be 

2*Be 

P3 

*Rp3 

Rp3* 

*Rp3 

Rp3* 

*Rp3 

T 

= .5* 

ta*E 

Rn 

b ta * 

eta* 

si on 

*£ * f 

+ 2*E 

s*En 

s*Er 

ts* 

ts* 

T 

t = Et 



De L rs=De L rs+Rts 

m*LQG(N*M) 

M-1/N) 

/N-2+1/M*2) 

n=1 ] OR (Version=0] THEN L1 



eta^PS 

Beta*R] 

a*E 

ta*[Rnn*E+Rn*En) 
ta*E r 



(En+E*P3*Rn/R] 

[Er+P3SE/R) 2 * En+E * tP3 " 1] * Rn/R] * Rn+E * Rnn)/R) 
*(Err+2*P3*Er/R+P3*(P3-1 ] *E/R*2) 

EXP[-Beta*R) 



[Rnn*E+Rn*En] 

Er 

= THEN L1 5 

1+E] 

*Efac 
*E f a c 

Enn*Efac+2*En*2 
Err*Efac+2*Er*2 



s*E 
nn 
ts*Er r 



s*En 

ts*Ei 

s*Er 



o 



5BQD 
5B02 
5604 
5606 
560B 
5610 
5612 
5614 
5616 
5618 
5620 
5622 
5624 
5626 
5628 
5630 
5632 
5634 
5636 
5638 
57Q0 
5710 
5720 
5730 
5740 
5750 
5760 
5770 
5780 
5790 
5800 
5810 
5820 
5830 
5B4Q 



SUB CubictCf [*] ,Rt(«] ] 

R=-Cf[0 

Q=-Cf (1 

p = — Cf r p j 

A=[3*Q-p-2)/3 

B=(2*p-3-9*P*Q+27*R)/27 

Fc1=A*3/27 

Fc = B~2/4+Fd 

IF Fc<0 THEN L89 

PRINT "Error-005" 

STOP 

L89: Psi=-B/t2*SQR[-Fe1 }] 

DEG 

Psi=ACS(Psi 1/3 

Fc=2*SQR{-A/3) 

Rt 1 =Fc*C0S(Ps1 J-P/3 

Rt 2)=Fc*C0S Psi+12D)-P/3 

Rtl3}=Fc*C0SlP8i+240)-P/3 

RAD 

SUBEND 

SUB Normod(H[*) .Ev [*) .Ft* 

DIM GM:2.1 :2),Gt(1:2 1:2 

FOR 1=1 TO 2 

Gh,IJ=-H(1 ,2)/[H[1 ,1 ]-Ev(I) ) 

NEXf I 

MAT Fg=F*G 

MAT G?=TRN(G) 

MAT Gfg=Gt*Fg 

FOR I=T TO 2 

Lt2 f IJ=SQRtEv(I)/GfgtI f I) ) 

LM ,li=G(1 ,I]*L(2,I7 

NEXf I 

MAT Li=INV(L) 

SUBEND 



,L1(*n 

»Fg(1 :2,1 :2) ,Gfg(1 :2,1 :2) ,LM :2,1 :2 



On 



6DDD 
60D5 
6010 
601 5 
6020 
6025 
6030 
6035 
6040 
6045 
6050 
6055 
6060 
6065 
6070 
6075 
60B0 
6085 
6090 
6095 
6100 
6105 
6110 
6115 
6120 
61 25 
6130 
6135 
6140 
61 45 
61 50 
61 55 
6160 
6165 
6170 
6175 
61B0 
61 85 
6190 
61 95 
6200 
6205 
6210 
6215 
6220 
6225 
6230 
6235 
6240 
6245 
6250 
6255 



! 6 point Gaussian Lengendre 



SUB TunLlVI ,V2 , F. M.T.Gam) 

OPTION BASE 1 ' ' ' 

DIM X[6),W(6) 

RdR A N = 4 3 TO 1 eP 8B083 »- 661 2093B6466,. 93 246951 4203 

READ X{N) 

X[7-N)=-X[N) 

NEXT N 

PAJA 467 91 393457 3 ,.360761 573048 ,.1 71 32449237 9 

rUn N — 4 lU b 

READ WIN] 

W(7-N)=W[N] 

NEXT N 

H=6.6234E-27 

K=1 .3B033E-16 

Kt=K*T 

A=V1-V2 

V1 h = SQR[V1 ) 

V2h=SQR(V2J 

B=[V1h+V2h) "2 

Pi2=2*PI 

L=Pi2*SQR(-2/F)/(1/V1h+1/V2h) 
C=H-2/(B*M*L-2J 

DEF FNCDSh(Z)=.5*(EXP(Z)+EXP(-Z) ) 

De L ta= [ B-C )/C 

IF DeLta<0 THEN L1 

Dfac=FNCosh(PI*SQR[DeLta) ) 
GOTO L11 

L10: Dfac=COS[PI*SQR(-DeLta)] 

LI/I: IF V2> = V1 THEN E0=-V1/Kt 

IF V1 >V2 THEN E0=-V2/Kt 

Vav=.5*(V1+V2) 

i&l = i C o lL0G(2 * M+Dfac)/ - Q14)/p i 2 ^2-Vav]/Kt 

CD c — d , d 

Eb=MIM(Eb1 ,Eb2) 

Em=.5*(Eb-E0) 

Ep=.5*(Eb+E0] 

Gam = 

FOR N=1 TO 6 

E=Em*X(N)+Ep 

Kt e= K t *E 

ALph1=PI*SQR( (Kte+VI )/C) 

ALph2=PI*SQR([Kte+V2)/C] 

Facp=FNCosh(ALph1+ALph2) 

Facm=FNCosh(ALph1-ALph2] 

Y=Ke*EXP[-E] 
Gam=Gom+W[N) *Y 
NEXT N 
L50: Gamfac=EXP(-Eb) 

Gam=Em*Gam + Gatnfac 
SUBEND 



6500 SUB PLtMOptI*) ,Nm f T(*] ,Lk I*] ,E[*) ) 

6502 DEG .„..,„„,. 

6504 PLOTTER IS 1 3 , "GRAPHICS" 

6506 GRAPHICS n nr , 

6508 LOCATE 19,119.18,95 

6510 DIM XtlcsM :11J ,YticsM:11 J 

6512 DATA . 01 , . 02 , . 025 , . 05 , .1 , .2 , .25 , . 5 , 1 ,2 , 2 .5 

6514 MAT READ Xtics 

6516 Xmin = 1QO0/T 1 ) , 

651 8 Xmax=1 0O0/T[Nm] 

6520 Xspan=Xmax-Xmin 

6522 FOfi 1=1 TO 11 

6524 Itab=I , # - % * , „ -r-t.r-t., inn 

6526 IF INTlXspan/XtlcsII J)<=6 THEN L30 

6528 NEXT I . . u , 

6530 L30: Xt i c = Xt i cs ( I tab , 

6532 Xstart=INT[Xm1n/Xticj*Xti c 

6535 IF'lSstSJixHSjf AND CXslop-Xma x >1 E-2*Xt i c ) THEN Xs top= ( I NT [Xmax/Xt i c ) +1 J *Xt i c 

6536 Ymin=Ymax=Lk[1 J 

6538 FOR N=2 TO Nm , i% 

6540 IF Ym1n>LMN| THEN Ym1n=Lk N 

6542 IF Ymax<Lk(N) THEN Ymox=Lk(N) 

6544 NEXT N 
S 6546 Yspan=Ymax-Ymin 
to 654B MAT Y t i cs=Xt i cs 

6550 FOR 1=1 TO 11 

6554 IF Q INT[Yspan/Ytics(I) )<=6 THEN L40 

6556 NEXT I 

6558 L40: Yt i c=Yt i cs ( I tab | t 

6560 Ystart=INT Ymin/Ytic *Ytic 

6 56 2 Y8top=(INT[Ymax/Yt1c5+1 )*Ytic 

6564 SCALE Xst a r t , Xstop , Ysta r t , Ystop 

6566 LINE TYPE 3 v «. v * «. 

6568 GRID Xt i c , Yt i ic ,Xst a rt , Ysta r t 

6570 LINE TYPE 1 

6572 AXES Xt i c f Yt i c f Xs ta r t , Yst a r t , 2 , 2 ,6 

6574 FRAME 

6576 CSIZE 3 

6578 LORG 8 «Ti-n \J±. • 

6580 FOR Ypos=Ystart TO Ystop STEP Ytic 

6582 MOVE Xstart,Ypos 

6584 LABEL USING "MDD . DDX" j Ypo a 

6586 NEXT Ypos 

6588 LORG 6 t „_„ v, . 

6590 FOR Xpos=Xstart+Xti c TO Xstop STEP Xtic 

6592. MOVE Xpo s , Ysta r t-Yt i c/1 

6594 LABEL USING "Z.DD«»;Xpos 

6596 NEXT Xpos 






6598 


LDIR 90 


66DD 


CSIZE 4. 


6602 


LORG 1 


6604 


ON 0pt(2 
L1 : ON 


6606 


6608 


La1 i MOV 


6610 


LABEL US 


6612 


fi #% "■* g-% ■ 

GOTO Lou 


6614 


La2; MOV 


6616 


LABEL US 


661 a 


GOTO Lou 


6620 


La3 : MOV 


6622 


LABEL US 


6624 


GOTO Lou 


6626 


La4: MOV 


6628 


LABEL US 


6630 


GOTO Lou 


6632 


L2: ON 


6634 


Lb1 : MOV 


6636 


LABEL US 


6638 


GOTO Lou 


6640 


Lb2: MOV 


6642 


LABEL US 


6644 


GOTO Lou 


6 6 46 


Lb3; MOV 


6648 


LABEL US 


6650 


GOTO Lou 


6652 


Lb4: MOV 


6654 


LABEL US 


6656 


Lout: LD 


6658 


MOVE Xst 


6660 


LABEL US 


6662 


MOVE T(1 


6664 


FOR N=1 


66 66 


PLOT 100 


666B 


NEXT N 


t? c? tj n 
ODD tJ 


LINE TYP 


6670 


MOVE 100 


6674 


FOR N=1 


6676 


PLOT 100 


6680 


NEXT N 


66B2 


RAD 


6684 


SUBEND 



) GOTO L1, 

pt(4) GOTO 

E Xstart-X 

ING "K";"L 

t 

E Xstart-X 

ING »K" ;»'L 

t 

E Xstart-X 

ING "K";"L 

t 

E Xstart-X 

ING "K»*;"L 

t 

pt(4) GOTO 

E Xstart-X 

ING »K" ;"L 

t 

E Xstart-X 

ING "K" ;"L 

t 

E Xstart-X 

ING »K" ;»L 

t 

E Xstart-X 

ING »K" ;"L 

IR 

art+ ,43*Xs 

ING »K" ;"1 

J.Lktl ] 

TO Nm 

O/TIN] , Lk( 

E 3 

O/TM ) ,Lk( 
TO Nm 
0/T(N) ,LM 



L. 8 I ■ La 

span* .1 
n [ k ) c 

span*.1 
n ( k ] c 

span* .1 
ntk] L 



span 
nlk] 



.1 
L 



Lb 1 f Lb 
span * .1 
og ( k) 



2 ,La3 ,La4 

5 f Ystart+Yspan 

c/mo La-s" 

5 , Ystart + Yspan 
c/mo Lacu L a-s" 

5 , Ysta rt + Yspan 
i ters/mo La-s" 

5 , Ysta rt + Yspan 
i ters/mo L ecu L e 

2,Lb3,Lb4 

5 f Ysta rt + Yspan 

cc/mo Le-s" 



span*. 1 5 ,Ysta rt+Yspan 
og(k) cc/mo Lacu La-s" 

span *.1 5, Ysta rt+Yspan 
og(k) L l ters/mo Le-s" 

span*.1 5,Ystart+Yspan 
og[k] L i ters/mo Lecu L 

Ban,Ystart-Yspan*.12 
00/T" 



N] 

1 ]-E(1 ] 
N]-E(N] 



*.27 

*.1B 

*.1B 

♦.12 
-s" 

* . CD 

*.16 
*.16 
*.08 

B-S" 



ON 
C/1 



BBOO 
68Q1 
B81D 
682D 
6830 
6831 
685D 
6860 
6870 
6880 
6890 
6900 
6910 
6920 
6930 
6940 
6950 
6960 
6970 
6980 
6990 
7000 
EOF 



SUB FlU0pt2,N f T(*],Xob(*] ,E(*) r P { * ) , E2 , 0pt7 ) 

OPTION BASE 1 

DIM AIN.3J f At [3 f N] f At a (3, 3 3 ,Atai [3,3) ,Atx(3) , 

FOR 1=1 TO N 

ACI t 1 3=1 

A{1.2)=0 

IF 6pt2=2 THEN L10 

IF 0pt7=3 THEN AII.2 = 

A{I,3)=-1/M .987*111)] 

G0T6 L11 



=LOG(T(I 



L10: IF 0pt7=3 THEN A ( I , 2 =LGT [T ( I 

All 3]=-. 45429/(1 ,9B7*Tlf J) 

L1 1 : N EXT I 

MAT At=TRN(A) 

MAT Ata=At*A, 

MAT Atai=INV(Ata ) 

MAT Atx=At*Xob 

MAT P=Atai *Atx 

MAT X=A*P 

MAT E=Xab-X 

E2=D0T(E,E) 

SUBEND 



/ 



